NGINEER: NG + 
uPRRY | JUL 1 4949 


THE QUARTERLY JOURNAL OF 


MECHANICS AND 
APPLIED 
MATHEMATICS 


VOLUME II PART 2 
JUNE 1949 


OXFORD 


AT THE CLARENDON PRESS 
1949 


Price 12s. 6d. net 





THE QUARTERLY JOURNAL OF MECHANICS) 
AND APPLIED MATHEMATICS i 


Editorial Board 


S. GOLDSTEIN R. V. SOUTHWELL 
G.I. TAYLOR G. TEMPLE 
together with 
H. JEFFREYS 
J. E. LENNARD-JONES 
N. F. MOTT 
W. G. PENNEY 
A. G. PUGSLEY 
L. ROSENHEAD 
ALEXANDER THOM 
A. H. WILSON 
J. R. WOMERSLEY 


Executive Editors 
G. C. McVITTIE Vv. C. A. FERRARO 


THE QUARTERLY JOURNAL OF MECHANICS AND APPLIED MATHEMATICS i§ © 
published at 12s. 6d. net for a single number with an annual subscription © 
(for four numbers) of 40s. post free. 


NOTICE TO CONTRIBUTORS 


1. Communication. Papers should be communicated to one or other of the Executive 5, 
Editors, by name, at King’s College, Strand, London, W.C. 2. 


2. Presentation. Manuscripts should preferably be typewritten, and each paper should ~ 
be preceded by a summary not exceeding 300 words in length. References to literature 
should be given in standard order, author, title of journal, volume number, date, 
These should be placed at the end of the paper and arranged according to the orderof 
reference in the paper. % 


3. Diagrams. The number of diagrams should be kept to the minimum consistent with 


clarity. The lines of the figures should be drawn in ink either on draughtsman’s paper 


or on good quality white paper. Each individual line in the figure should bear reducing ~ 
to one-half of the size of the original, and great care should be exercised to see that the 2 
lines are regular in thickness, omg where they meet. Lettering of the figure should ~ 
be in pencil and should be sufficient to define clearly the lines and curves in it. The 
writing of formulae or of explanations on the diagram itself should be avoided. All 
explanations of symbols, etc., should be given in underline. Contributors should indi- 
cate on their manuscripts where figures should be inserted. 


4. Tables. Tables should preferably be arranged so that they can be printed with the 
columns parallel to the longer edge of the page. 


5. Vector Notation. All single letters used to denote vectors in the manuscript should 
be marked by underlining with a wavy line. Scalar and vector products should be 
denoted by a. and a a b respectively. 


6. Offprints. Authors of papers will be entitled to 25 free offprints. 
7. All correspondence other than that dealing with contributions should be addressed 
to the Publisher: 
GEOFFREY CUMBERLEGE 
OXFORD UNIVERSITY PRESS 
AMEN HOUSE, LONDON, E.C, 4 








form 


and t 





ON A FAMILY OF ROTATIONAL SPATIAL GAS FLOWS 
By R. PRIM and P. NEMENYI 
(Naval Ordnance Laboratory, White Oak, Maryland, U.S.A.) 
[Received 18 May 1948] 


SUMMARY 

An investigation is made of steady flows of an ideal gas in the absence of body 
forces for which the three velocity components in a cylindrical coordinate system 
based on a plane isometric net depend on only one isometric coordinate. It is shown 
that such flows exist only if the basic isometric net consists of logarithmic spirals 
or their limiting cases (concentric circles, radial straight lines, and parallel straight 
lines). Several large classes of rotational and spatial gas-flow solutions are obtained, 
some in explicit forms involving an arbitrary function, and others as solutions to 
ordinary differential equations involving several arbitrary parameters. 


Introduction 

ToriMIEN (1) has investigated plane, irrotational, steady flows of an ideal 
(ie. thermodynamically perfect, non-viscous, and thermally non-conduct- 
ing) gas in the absence of body forces for which the velocity components 
with respect to an orthogonal coordinate net are functions of one coordinate 
only. In the present paper these conditions are generalized by the addi- 
tion of a velocity component normal to the basic coordinate net and by 
the removal of the restriction to irrotational flows. At the same time, 
however, the conditions are specialized by the consideration only of those 
plane orthogonal nets which are isometric, that is, which could be obtained 
by conformal mapping from a Cartesian net. 


Formulation of the problem 

We employ the formulation by Munk and Prim (2) of the basic equations 
governing steady flow of an ideal gas in absence of body forces. This 
formulation yields the continuity equation 

div| (1— W2)¥7-DW] = 0 
and the integrability condition for the dynamic equation 
W /\ curl W 
a 

where y denotes the adiabatic exponent and W the reduced velocity vector 
defined in terms of the actual velocity vector V and the ultimate velocity 
Magnitude a (constant along any individual streamline) by 


= 0, 


curl | 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 2 (1949)] 
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For the present investigation, (1) and (2) will be referred to a cylindrical 
coordinate system based on a plane isometric net. The squared element 
of are length in this system is given by 

ds = g°(&, n)(dé?+-dy*) +dz*, (4) 
where the limitation of the (€,) net to a plane requires that logg be a 
harmonic function of € and », that is, that 
67 log 67 Ic ‘ 
ete! = Velogg = 0. (5) 
0&? or” 


teduced velocity fields of the following type will be considered: 
W= EB, u(€) t M1 V(E)+2Z, 9(€), (6) 
where &), y, and Z, are unit vectors in the €, y, and z directions, respec- 
tively. 
Discussion of the equations 


Specialization of (1) and (2) by (4) and (6) yields four equations (in 
addition to (5) ) restricting the functions g(€, 7), u(€), v(€), and q(€): 











v = (log g)+-u A log| gu(1— W2)¥- Wi == @, (7) 
on 0€ ; 
é | u(dq 6) _ 6 [ u(dg/dé)] _ (8) 
e&|g(1—W*) | énlg’—W*)} | 
and v ‘ (logQ)+u Ss log ; — — 0, (9) 
where Q = a + os (log g)—u (log g). 


The double equation (8) requires that one or more of the conditions 
u = 0, dq/dé = 0, or ég/én = 0 be satisfied. When the various possible 
combinations of these conditions are examined in the light of the restric- 
tions imposed by (5), (7), and (9), four distinct cases emerge. These cases 
are characterized by the conditions 
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It should be pointed out that, to be physically meaningful, the magnitude 
f the reduced velocity vector must be between zero and one. Choice of 
‘arbitrary’ functions and constants in the following analysis must always 
he made consistent with this limitation.) 


Case I: wu = 0, eg/én #0 

For u 0, 6g/en ~ 0; (7) requires that v be also zero. No restriction 
s then imposed on q(&) or g(&,) by (9), and we are left with a uni- 
lirectional flow along the z-axis with the velocity magnitude an arbitrary 
function of € where the (€,) net may be any isometric net. This result 
is relatively trivial, being a special case of the well-known family of flows 
parallel to the z-axis, with W an arbitrary function of x and y. 


Case II: wu = 0, ég/en = 0 

For u = 0, ég/én 0; (7) and (9) impose no limitation on the functions 
é) and g(é). However, (5) limits the permissible isometric nets to those 
or which logg = A€+C, where A and C are constants. For A + 0, the 
lines € = constant form a set of concentric circles; for A = 0, a set of 
parallel straight lines. 

Thus we have, basically, flows whose streamlines are helices lying on 
concentric circular cylinders and in which the tangential and axial velocity 
“mponents are arbitrary functions of the radius. In the limiting case, 
4= 0), the cylindrical surfaces are a family of parallel planes in each of 
which the velocity vector may be chosen arbitrarily. 

Case III: wu 4 0, dg/dé = 0 

For this case, the streamlines no longer lie in the surfaces € = constant, 
but the z-component of velocity is a constant, say q(é) = Q. 

If L denotes the function — (log g) it can be shown by differentiation of 

Cc 7) 

)), (7), and (9) with respect to » that the following equations are necessary 

restrictions of the function: 


C2 CL 2 CL 
(log (I y— 0, 10 
sen =e) 3 Ox? one) (10) 
OL o? ( 
—(I og )=; < (lee) = (11) 
C&C 7 \ on 07) 
ae y 0. (12) 
og* o7n* 


this set of necessary conditions, together with (5), leaves the following 
possibilities for the function log g: 


logg = cP cos F(é—D)+AE+ Bn tO (F #0) (13) 
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or, logg = $G(n?—&)+ Hén+AE+ Bn+C. (14) 


However, these conditions have been obtained by using (7) and (9) only 
in their differentiated forms and hence represent necessary, but not 
sufficient, limitations on logg. The undifferentiated forms of (7) and 
(9) further restrict (13) and (14) to 


logg = AE+ Bn+C. (15) 


Hence, the lines € = constant of the basic coordinate net form a set of 
logarithmic spirals or their limiting cases (A = 0, radial straight lines: 
B = 0, concentric circles; A = B = 0, parallel straight lines), 

Under this restriction of logg and of q, (7) and (9) reduce to a pair of 
ordinary differential equations restricting u(&) and v(é), namely, 


Bv+Au- Fe loslu( -Q?—u?—v?)"y-] = 0 (16) 
lv ) 

and ay atl mal re we. (1—Q?—u?—v?)., (17) 
dé a )) 


These equations determine a five-parameter family of spatial gas-flow 
solutions, only very special cases of which were known before. For D = 0, 
Q ), the reduced velocity field is irrotational and plane, and Tollmien’s 
flows ce on a logarithmic spiral coordinate net (including the Bateman- 
Taylor spiral flows as special cases) are obtained. For A = 0, D= 0, 
Poritsky’s (3) generalization of the irrotational Prandtl—Meyer corner flow 
is obtained. For A = 0, D + 0, a generalization of the Prandtl—Meyer 
corner flow to rotational flows is obtained. A detailed study of the latter 
case in an independent investigation by Prim (4,5) has yielded, among 
other results, a simple family of explicit solutions for rotational gas flows. 


Case IV: u + 0, ég/én = 0 


As in ease II, the lines € = constant of the basic coordinate net are 
restricted to families of concentric circles or of parallel straight lines. 
However, in the present case the streamlines do not lie in the surfaces 
€ = constant and the z-component of velocity is not necessarily a constant. 

For A + 0, that is, the surfaces € = constant forming a set of concentric 
circular cylinders, the constant C may be taken as zero without loss 
of generality, and the radius of the cylinders € = constant may be intro- 
duced as independent variable by the substitution 


r = eA€, (18) 


The restriction of the velocity components u(r), u(r), and q(r) by (7), (8 
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and (9) can then be satisfied by the introduction of the auxiliary dependent 
variable U(r) through the definition 


r 


U(r) ( p| | u*(p) v*(p) q?(p) |v” 1) dp (19) 


r 


if the velocity components are given in terms of the function U by 


y= tri on 
J 
a rit 2] 
] (21) 
K 
ind ( U+L 22 
] ] (22) 


capital letters represent arbitrary constants). Compatibility of the 
equations (19) and (20)-(22) requires that U be a solution of the ordinary 


differentia! equat ion 


- (av 1)/) J2p-2 v9( 2) - vr 20244 u+L) = 1, 


di dr 
(23) 


Equations (20)—-(23) determine a five-parameter family of rotational spatial 
gas flows. A common property of these flows is a helicoidal symmetry— 
any streamline of a given flow can be brought into coincidence with any 
other by a translation along the z-axis and a rotation about it. It should 
be added that by going back directly to the formulae (7, 8, 9) we can 
reduce the problem to a differential equation of the form du/dr = ¢(r, u), 
but the function ¢ is rather complicated. 

For A = 0, the surfaces € = constant form a family of parallel planes 
parallel to the z-axis. Without loss of generality we may refer our 
lows to a Cartesian coordinate system by taking £ = xand y = y. Again 
tis convenient to handle the limitations imposed on the functions u(7), 

v), and g(a”) by (7) and (9) through the introduction of a new dependent 
ariable by the definition 


V(x) ( l u=(s) v(s) q?(s) |v” ) ds. (24) 


The velocity components are then given in terms of V by 


1V\-U 
u — MI ¢ 
FS, \ 


to 
or 
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IN os 
v U if (26) 
P 
and ( V+_R. 97 
un f] i (21) 


Compatibility of (24) and (25)-(27) then requires V to be a solution of the 
ordinary differential equation 


dV\(y-Diy dV\-2” N? P \2 
. + M 45 yay LF ‘ 
(Te) bs 4 M2 ( V ] ] (28) 


a quadrature, if we go back directly to equations (7), (8), and (9). Under 
the present assumptions these equations take on the form 


u = A,(1—W?2)-"7-», (7! 
uq’ = C,\(1—W?), (8') 
uv’ = B,(1—W?). (9) 


From (8’) and (9’) it follows that 


Combining this with (7’) we obtain v in terms of w: 


C, D, [ODsB 


- ‘i i , Bi(i—Dj—u*— Az tu-))} 
By{1+(C3/ BY} ~~ L(Ci+ Bi? B+ 03 


From (9’) we obtain now: 


dv ‘ B, AY 1 1) 
U u >, AS or, 
du iol 
, dau B, A*-1u-v 
and hence u <— # ion 3 
dx dv/du 
: i dv 
therefore a= Bi Ai’ | «7 i du. 
du 


Integration by parts yields finally the relation 
0 By-2 
B, Al-!2 = wvv+ ; - a se 
(2y—1)(3y— 2)(B2-+ C2) 
7 eval CiDi Bi, BY I—Dj—wW—Ay ‘u-9-) 18 
(BEC) B+} | 


ny 


In this particular case it is possible, alternatively, to obtain w in terms of 
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ON A FAMILY OF ROTATIONAL SPATIAL GAS FLOWS 135 
Remarks on the results of the investigation 


A major result of the above investigation is the manner in which the 
choice of the basic isometric net is limited by the governing gas dynamic 
relations. The above analysis establishes the following: 

THEOREM. T'he only plane isometric nets for which it is possible to find 
steady flows of an ideal gas in the absence of body forces, such that the velocity 
components with respect to cylindrical coordinates based on the net are 
functions of only one isometric coordinate, are families of logarithmic spirals 
and their limiting cases. 

This result has been established earlier by the present writers in a 
preliminary publication on this and reiated questions of gas dynamics (6). 
It should be noted that, as Tollmien’s analysis (1) has shown, exactly the 
same isometric nets are admissible as coordinate systems for the special 
ease of irrotational flow. 

Apart from the general theorem, and from the previously known flow 
solutions occurring in cases I and IT, the analysis has yielded explicitly or 
in the forms of solutions to ordinary first-order differential equations 
several broad families of new gas flow solutions, including truly spatial 
rotational flows. The detailed study of representative examples of these 
slutions can be expected to yield valuable insight into the nature of 
rotational gas flow. 
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SUPERSONIC FLOW PAST THIN WINGS 
‘ % T al _ 7 
I. GENERAL THEORY 
By G. N. WARD 
(Department of Mathematics, The University, Manchester) 
[Received 10 August 1948] 
SUMMARY 
The linearized potential problem of supersonic flow is solved for nearly plane 
wings of arbitrary shape of planform and section. The velocity field is divided into 
its symmetrical and antisymmetrical parts, which are treated separately. The 
results for the symmetrical part have been given before, by A. E. Puckett, and are 
only considered briefly. The results for the antisymmetrical part are determined by 
extending an integral equation method used previously by J. C. Evvard in a more 
restricted problem. The Kutta—Joukowski condition is applied to the velocity at 
the trailing edge in order to make the solution determinate. Expressions are given 
for the potential on the wing, from which the velocities can be obtained. For the 
smaller aspect ratios, the analysis becomes complicated and only an indication of 
how the results can be obtained is given. The velocity is singular at the subsonic 
leading edges and the forces on the wing have to be determined by a limiting process 
which is illustrated by ‘obtaining expressions for the lift and drag forces. All th: 


results are completely general and no special cases are considered. 


1. Introduction 

THE problem of determining the supersonic flow of a perfect gas past a 
thin wing set at a small incidence to an initially uniform stream is solved 
by assuming that the flow is inviscid and isentropic, and that the perturba- 
tion velocities caused by the presence of the wing are so small that their 
squares and higher powers may be neglected in the equations of motion 
The wing surfaces are assumed to be nearly plane, so that wings with a large 
dihedral angle are excluded. The problem is solved by dividing the velocity 
field into its symmetrical and antisymmetrical parts and treating each 
independently; the results are then superposed to give the complete solu- 
tion for the general case. The symmetrical part corresponds to a wing 0! 
small finite thickness set at zero incidence which experiences no lifting 
force: the solution in this case presents no difficulty and is obtained 
immediately from the general solution of the equation of motion. The 
antisymmetrical part corresponds to a wing of infinitesimal thickness and 
the flow is determined by using an integral equation method following ® 
transformation to characteristic coordinates. It is assumed that the flow 
over the trailing edge of the wing is smooth in accordance with the 
hypothesis of Kutta and Joukowski, and a trailing vortex sheet is formed 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 2 (1949)] 
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SUPERSONIC FLOW PAST THIN WINGS 137 
behind the wing: without this assumption the solution of the anti- 
symmetrical problem is indeterminate. The velocity components exhibit 
singularities at the subsonic leading edges and the aerodynamic force and 
moment on the wing have to be obtained by using a somewhat elaborate 
limiting process, which is illustrated by the determination of the lift and 
drag components of the force. 

The solution for the symmetrical case has been given previously by 
Puckett (ref. 1), and, recently, Evvard (ref. 2) has solved the problem 
of obtaining the flow pattern near the wing surface for a restricted class 
of wing tips with curved subsonic leading edges. The use of the trans- 
formation to characteristic coordinates is due to Evvard, and the present 
treatment extends this idea to solve the completely general problem. 


2. Mathematical formulation 

Let x,y,z be a right-handed system of Cartesian coordinates such that 
r is measured in the direction of the undisturbed stream, which has a 
velocity U and Mach number & relative to the coordinate axes, and such 
that z is measured vertically upwards. Then, to a linear approximation, 
the components w,v,w of the perturbation velocity are the components 
of U grad d, where ¢ satisfies 


> 7 
m 24 C-@O c 2h 


Be”. Ca mee O, B M?—1., (1) 
ox oy~ cz” 

The perturbation pressure, p, is given, to the same order of approxima- 
tion, by the linearized Bernoulli equation 

p pU? od Ox, (2) 
» being the density of the undisturbed stream. 

The general solution of eqn. (1) for z > 0, say, when the boundary data 
are specified on z 0, as they will be in our case, can be obtained from 
Hadamard’s solution of the general second-order hyperbolic differential 
equation (ref. 3) by using the reflection method (ef. ref. 3, § 156) and may 
be written in the alternative forms 











ee | (op(x", y’,2’)) dx'dy’ as (3) 
a Jy | 2’ Jona9n (e—a’)?— B2(y—y’)*— B*2?} 
O! 
ere 
Or. y, ) A( TF OS 0) : ‘ <_< 719 4229.9) dx'dy’, 
7 J. dz|/{(a—a’)?— B*(y—y’)*— Bz*} pili 
(4) 


where, in both forms, the integration is taken over the area of the plane 
0 for which 

f ; > | V2 1 52 r 

’ a— BY (y—y’)? +27} (5) 
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and the star (*) in eqn. (4) denotes that the ‘finite part’ of the improper 
integral is taken (cf. ref. 3). The form of the inequality (5) is an expression 
of the fact that disturbances are not propagated upstream in supersonic 
flow. 

We restrict ¢ to be continuous except through the wing and the trailing 
vortex sheet. Sufficient conditions for the validity of eqns. (3) and (4) are 
that (@¢/éz),_, should be integrable over any finite area of the plane x = () 
and should possess discontinuities and singularities on only a finite number 
of ares in the plane, and that 4(~, y, +-0) should be continuous and bounded. 
These conditions are not necessary but are adequate for our purpose. 

Consider a wing whose surfaces lie approximately in the plane z = 0. 
Let the vertical projections of the wing and its trailing vortex sheet on 
the plane z = 0 be denoted by S and T respectively, and let the remainder 
of the plane z = 0 be denoted by R. Let 1,,m,,n, and /,,m,,n, be the 
direction cosines of the outward normals to the upper and lower wing 
surfaces respectively, so that 1,,m,, /,,m, are small and »,,, are approxi- 
mately equal to unity. 

The boundary condition of zero normal velocity at the upper wing 
surface is 


l 


(U+u)+m,v+nyw 0. (6) 

To our degree of approximation, the boundary condition may be applied 
on the plane z = 0, instead of on the actual wing surface, and similarly 
for conditions on the vortex sheet. Then, by neglecting small quantities 
of the second and higher orders, the boundary condition (6) becomes 


w od os 
ii = (=) —l,, (4) 
C% }2=+0 


and similarly, for the lower wing surface, the boundary condition is 


= == (-*) l (8) 
l 02 | -=-9 


where z +0 and z 0 denote the upper and lower sides of the 
plane z = 0 respectively. 

If 6 = d,+¢,, where ¢, is an even and ¢, an odd function of z, then 
dy is the potential for the symmetrical part of the velocity field and 4, is 
the potential for the antisymmetrical part, and the boundary conditions 
for ¢d, and ¢, over the area S are 


(“=*) - (2) = —}(1,+1,) (9) 
G2 Je, 62 }--— oO 


and (=) (=) — (10) 
02 Jeni 02 | -n~9 


respectively. 
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For the symmetrical problem, @¢)/éz is an odd function of z, and, since 
the z-component of the velocity must be continuous everywhere except 
through the wing, we have the following boundary value problem for ¢o 
satisfying eqn. (1) in the region z 0, say: 

(i) éd,/ez takes given values on S$ (from eqn. (9)) ; 
(ii) Od,/ez 0on Rand T; (11) 


0 

(iii) ¢, is continuous everywhere in z > +0. 

The potential dy is given everywhere in z > 0 by applying eqn. (3), 
since @¢,/€z is known everywhere on z = 0, and this solution satisfying 
the conditions (11) is known to be unique. 

The potential in z < 0 is also determined since ¢, is even in z. This is 
the result given by Puckett (ref. 1). 

For the antisymmetrical problem, 0¢,/éz is known only on S. The 
potential d, and é¢,/éx and 0¢,/éy are odd in z, and since the pressure 
must be continuous everywhere except through the wing, we must have 
é¢,/ex = 0on Rand T. Furthermore, ¢, is continuous through R so that 
6, =0o0n R. In addition it is assumed that the Kutta—Joukowski con- 
dition is satisfied at the trailing edge: this condition is that velocity 
remains finite at the trailing edge. It is shown later that this condition 
implies that grad, is continuous at subsonic trailing edges, but it may 
be discontinuous at supersonic trailing edges. Thus the boundary value 
problem for ¢, satisfying eqn. (3) in z > 0, say, is: 

(i) €¢,/éz takes given values on S$ (from eqn. (10)) ; 

(ii) 0¢,/éa OonT: 

; QonR: (12) 


iv) |gradd,| is finite on the common boundary of S$ and T; 


(iii) db 


(v) ¢, is continuous everywhere in z > +0. 


It is assumed that the solution for ¢, under these conditions is unique. 

The potential 4, in z < 0 is also determined since ¢, is odd in z. 

If the subsonic trailing edge of the wing is not sharp and the Kutta— 
Joukowski condition cannot be applied, then the solution for ¢, is not 





unique, and the flow pattern is not determinate on the present assumptions. 
It is sufficient to determine the potential ¢, on S, since, from the con- 
ditions (12), 4, is then known on R and T, and is therefore determined 
everywhere by applying eqn. (4). We shall return to this point later, in § 5. 
In order to calculate 4, on S, we take the special case of eqn. (3) when 


0 to obtain 


| l ( f (ed dx'dy’ 
| b,(x,y, +0) | ' — (8) 
) 7 | | 02’ Jrasigy((2—2' P— By—y’)*} 
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where the integration is taken over the area of the plane z = 0 for which 
xv <*#tBy-y’). (14) 
We transform this integral by putting 
x = x— By, B = «+ By, (15) 
so that «,8 are characteristic coordinates in the plane z = 0. 
Now OXY) 1 (16) 
@(a,B) 2B 
; 1 [é¢ 
and by putting ~ = —f(a, B) (17) 
27 B\0z }.- +9 
B’) da’dp’ 
eqn. (13) becomes : [a 3P )dadp 18 
q ( ) beec es dy, [ iT v')(B— -p’) y? (13) 


the integration being taken over the area for which 


ye ee p’ < B. (19) 


The function f(a, 8) is not known on R and T, and our procedure is to 





Fic. 1. 


set up integral equations to determine it there. In many cases it is not 
necessary to solve these equations completely in order to determine 4, 
on S, as was found by Evvard. 


Let C be the boundary of the region S and, in Fig. 1, let L, M,N, P be the 
points on C at which the tangents LL’, MM’, NR, PR are characteristic 
lines for eqn. (1) in the plane z = 0, making the Mach angle with the 


a-axis. LM is then the supersonic leading edge of S. LQ, MQ, NN’, PP 


are % 


tang 
vort 
Le 
let 


Let 
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are also characteristic lines. 7, S,, 7, S, are lines parallel to the x-axis, 
tangent to C at 7, 7, and are the boundaries of T, the projection of the 
vortex sheet. 

Let O, the intersection of LL’, MM’, be the origin of coordinates, and 


let L be the point (ap, 0), 
M be the point (0, Bp), (20) 
T,, be the point (a, 8;), 
T,, be the point (ag, B.). 
Let the equation of C be specified as follows: 
PLM: B = A,(a); MNP: B = A,(a) (21) 
or LMN: «= B,(8); NPL: «= B,(B) (22) 


so that A,, As, B,, B, are single-valued functions on the ares of C for which 


they are defined, if, for simplicity, we assume that the curvature of C has 
the same sign everywhere, as depicted in Fig. 1. If this condition is not 
satisfied, the following analysis may require modification, but this is 
easily carried out. 

The potential in the region LQM of S may be calculated immediately. 
From condition (12, iii), an application of eqn. (4) shows that ¢ and all 
its derivatives vanish upstream from the envelope surface of all the 
lownstream characteristic cones having vertices on the supersonic leading 
edge LM, which we call the leading characteristic surface. Thus 

(€¢/€z).-9 0 
on that part of R which lies upstream from L’L.MM’, and an application 
of eqns. (17) and (18) determines ¢ on LQM since (0¢/éz),_ ,. is known over 
the relevant part of z 0. More generally, an application of eqn. (3) 
determines ¢ in the space between the leading characteristic surface and 
the downstream characteristic cones from L and M. 
We now proceed to calculate 4, over the remainder of S: this is most 


conveniently done in successive stages. 


3. Independent wing-tip flows 

We consider first the case when the aspect ratio is sufficiently large for 
the flows over the two wing tips to be independent, i.e. such that the 
intersection Q of the characteristic lines LQ, MQ in Fig. 1 lies outside S, 
and calculate d, on that part of § downstream from MQ. 

For the wing tip shown in Fig. 2 (which is lettered to correspond with 
Fig. 1) on S, let 


f(a, B) A(a, B); (23) 
In the region M MT, S, on R. let 
f(«, B) (x, B), (24) 
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and in the region N’T, N on T, let 


f(a, B) = p(x, B). (25) 


p and p, are, as yet, unknown functions. 





For a point («,f) in the region M’MT7, T, (i.e. B > A,(a); O0<a< 
we have from eqn. (18) and condition (12, iii) 





$y 





This equation is satisfied if 


»(ax) 


B Ayla 
f mind “PMB a _ 
J (8—P’) V(B—B’) 


A(x) A;(a) 


27) 


This is an integral equation of Abel’s type for (a«,f) (so also, of 


course, is eqn. (26) for the contents of the brackets) and its solution which 
we shall require later is (cf. ref. 4) 


Ag(a) 
rox, B’)/{A (x) —B"} ap’ , 
: siaasiipnilaanaa Sea nln A St : (28) 
Mah) = — Sg. AA)} B—B’ 


For a point («,8) on S in the region 7} 7, MQ (i.e. B,(8) < « < B,(B): 
By < B < B,), we have from egn. (18) 


B 
[ __ oe ges 
J v(a—a') \(B—B’) 
Bi(B) Axa’) 
B,(B) 


$y 


B A,(a’) : 
f= fue, B a | Ma’,B') dB'| (ag 
i. ———— ! — 

’ Vla—a')|_ 


~V(B—B’) J 


19(a’) A;(a’) 








Ki 
on $ 


or 
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From eqn. (27) the last term vanishes and we have finally, for points 


on Sin the region under consideration, 


Y B 
F da’ i A(a’, B’) dB 
db ; —- — (30) 
. V\o 7 | (B B’) 
“lad 
or (X,Y, 0) ! (a) dx ‘dy’ es (31) 
1 \02" J vaso {(e—2’)? Biy—y )*} 


‘ 


the integration in eqns. (30) and (31) being taken over the area ©, bounded 


by characteristic lines, shown in Fig. 2. This result is due to Evvard 
ref. 2), who obtained it by essentially the same method. 
It remains now to calculate the potential in the region NT, 7. 
On T, for z +-0, the potential is a function of the coordinate y alone 
ind we can therefore express it as 
¢, = F(a—8), (32) 
Ff being an unknown function. Since the potential is continuous over 


0, it will be zero on the line 7, S, and we have 


F,(a,—B,) = 9, (33) 
8, being the coordinates of 7. The line 7; S, has the equation 
x—B %,—B; = y,, say. (34) 


For a point on T in the region V7, NV, we can apply the result of eqn. (30) 
{T is considered to be part of the area for which A is known, and A is 


replaced by up, over this area. We have then, on T, 


B Ada’) 
> da’ ( [ p,(a’,B’) dp’ Aa’, B’) dB’| 
0, -: — F\(a—B). 
a (x x )] | a) | 3 B’) V(B B’) } if ; B) 
Axa Aj(a’) (35) 


This equation is an integral equation of Abel’s type for the contents 
the brackets, and the solution is 


} Lo x 
* u(x, B’) dp’ [ ii B’ ) dp’ | [ Fy(a' -B) da’ 
(B—B’) 





; ~_ 
; cK (36) 
d (B—B’) 7. 
Ao(a 43 (a B+y1 


where F’, is the first derivative of F,; and since F,(y,) = 0, eqn. (36) is 


s(a—a’) 


nother Abel integral equation for »,, and its solution is 


1o(a) 


l [ eee x )—B'} 
tmr/{B—A,(a)} | B—p’ 


4;(a) 


u.(a, B) 


ew 


l 


c dp’ { Fy(«'—B’) d 37 
7 Op , V(B—B’) oa — 


A(x) "+1 
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The unknown function F, has to be determined by the condition that 
pu, Must be finite at the trailing edge 7,.N, where B = A,(«). It will be 
seen that the second term in eqn. (37) must contain a term which will 
cancel the singularity at the trailing edge contained in the first term, 
We need not actually solve for F, in order to determine the potential on 
the wing. Let a function G,(&) be defined by 

g 


l * Fi(E') dé’ 
Gl) = - | e+ .. (3 
mms ME=€) 
v1 
Then the second term in eqn. (37) is 
B . 
12 f Gy(o—B') dp’ _ AP Gyfo— Ayla} fF Glo B') AB] og 
7B J y(B-B) = ie A) G- Bp’) fw 
Axa A2(a) 


As the trailing edge is approached, 8 - A,(«) and the first term in the 
right-hand side of eqn. (39) becomes infinite, while the second term 
vanishes, as may easily be verified on taking account of the continuity 
of the function F,. Thus, by combining equations (37) and (39), if p, is 











to be finite at 8B = A,(«) we must have 
Ax(a) " F 
G{a—A,(a)} = | aah a - (40) 
ii * . 
whence 
] a 1 , ] F ve Qt 1p’ 
p,(a,8) = —/{B—Ag(a)} rid icra GaP 
al J (B- Aja)—ps 7, v(B—B’) 
A;(a) A2(x) 
(41) 
We now calculate the value of pu, at the trailing edge. Let 
Ap(x) 
i 
J = \{p—A,(a)} [ (Ma, P )—Mo, Asay (42) 
7 (B -p’ WtA(a)—B} - 
A(x) 


and let the function A be such that 





lim J = 0. (43) 
Then, since —— 
: a As x ) 
V{B—A,(a)} | ox Wa Fi si 7—2tan-* || F or 
— (44) 
and from eqns. (41), (42), and (43), 
lim p,(a,B) = Afa, A,(a)}. (45) 


B—A2(«) 
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Equation (43) is certainly true if A satisfies the conditions 
(i) |A(a, B)—Afa, Ag(a)}| < K{A,(a)—p}® for A,(a) >B >k > A,(a), 
k 

(ii) | JA x, B)—A{a, Ag(x)} dpB< L, (46) 
Ax( x) 

where K, L, 5 are positive constants, and these conditions are more general 

than those usually required in any physical problem. 

This shows that @¢,/éz is continuous at the subsonic trailing edge, and 

hence, from eqns. (13) or (18), 4,, é¢,/éa, and 0¢,/0y are continuous there. 

Thus grad 4, is continuous at the subsonic trailing edge. 

The potential at a point (a,8) on S in the region NT, 7; (Fig. 2) is 

obtained by using eqn. (30), as was done when deriving eqn. (35), and is 


B 
F ; da’ r A(a’, B') dBr 
x’) | ; 








5 ~ A 
a J Ie J <(B—B’ 
BiB 1;(x’) 
B,(B) ; B ’ ’ Ax(a’) ’ , 
7 l da’ | [ mela’ ) dp’ | Aa’, B)aB'\ (47) 
J Vea) J B-B) V(B—B’) 


By using the result of eqn. (36) to replace the contents of the brackets 
in the second term of eqn. (47), and remembering the definition of G, 


eqn. (38)), we have finally 


x , Pree ae ee ; 
$b; Sage § See [ ———_—, 
J Va’) J BBY J 7 aa’) 
By(B Ay(a B+y1 


where G, is given by eqn. (40) in terms of A. This result may be compared 
with eqn. (30). 

On the subsonic trailing edge 7, the first term in eqn. (48) vanishes, 
and by substituting for G, from eqn. (38) and inverting the order of 
integration, it may be verified that the potential is continuous there. By 
lifferentiating eqn. (48) and using eqn. (40) it may also be verified that 


6,/e2 = 0 on the subsonic trailing edge. 


We have now found the potential all over the wing tip on z = +0, and 
the results will be the same but with reversed signs on z —0. The 


| Potential on the other wing tip may be written down immediately by 


ising the formal symmetry of the problem with respect to « and 8B. Thus 
the potential is known at all points on both sides of the wing and the 





elocity can be found by differentiation. The calculation of the potential 
ind the velocity for points not on the wing will be considered later. 

4. Independent subsonic edges 

, When the wing-tip flows are not independent in the sense of the last 


section, an application of the results obtained already enables us to calculate 
5092.6 


L 
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the potential on the wing for the case when the characteristic lines LQ, 
MQ intersect on S and intersect C in the supersonic trailing edge PN 
so that the flows at the subsonic edges are independent. Two possibilities 
arise: (i) when the downstream characteristic lines from 7, and 7, intersect 
outside S, and (ii) when these lines intersect inside $. We shall consider 
the latter case for the moment, since the results for the former may be 
obtained by omitting certain of the terms which arise from the presence 
of the vortex sheet. 
Consider the wing planform shown in Fig. 3 (a). 





Fic. 3 (a). Fic. 3 (b). 


From the results of the previous section the potential is known on $ 
except in the region DGQ. We shall calculate first the potential in the 
, FH, G, 
for the left-hand tip correspond to the functions p, p,, F,, G, respectively 
for the right-hand tip, and let y, = 8,—a«, so that the equation of the 
line 7, S, is B—a = yo. 

The potential for a point («,8) in the region E FQ’ can be written down 
immediately from eqn. (48), by substituting v and v, for A in the appro- 
priate regions, in the form 


region EFQ’ denoted by ‘a’ in Fig. 3(a). Let the functions v, v 


8 





S P Bp —* A;(a) ‘ Bi(B) ; 

d= | ~. } See j Me | G,(a'—B) dor , 

tS (a—a’)| J (B—B’) WB—B) |* J 7 e—o') 
B,(B) A;(x) A(x’) Boy 
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fda’ Fy,(a’,B") dB’ ” F" v(a’,B’) dB’) 
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On inverting the order of integration in the second term containing A 
and in the terms containing v and y,, and by using the results for the left 
tip corresponding to eqns. (27) and (36), we can eliminate v and py, to obtain 


x ; B B,(p) Aj(a) 
* da | A(a B’) dp’ [ dx’ [ A x’, B’) dp’ 
ra J alla v) | (p B’) 3 (x x’) ,(B—B’) 
BB Ay(x By{ Ay(a)} A;(a’) 
By(p) A;(a) 
[ G,(a’—B) da’ | & G,(B’—a) dp’ (50) 
+ _ ») 
: \(a— a’) ~ A(B —P’) 
B.(B) ; 
ij ~ ’,B)d 7 
where G,{8—B(p)} = | = a... (51) 
: : p | B,(B) “- 
B,(B) 7 
The first two terms of eqn. (50) may be written in the form 
| [ X(a’, B’ deani J f | A(a’, B’) da’dp’ (52) 
. 52 
Ji(a—a )(B ,% Vi(a—a’)(B B’)} 


er by characteristic lines, are 


>> 


where the areas of integration &,, &. 
shown shaded in Fig. 3 (b). The area &, may not exist in some cases of course. 

For the region ‘b’ in Fig. 3(@) we have to omit the integral containing 
@, in eqn. (50), for region ‘c’ we have to omit the integral containing G4, 
md for region ‘d’ we have to omit both of these integrals. The potential 
for region ‘d’ is given by the expression (52), and this result for such 
regions was found by Evvard (ref. 2). 


The general supersonic wing 

When the flows over the subsonic edges are not independent, the 
calculation of the potential everywhere on the wing is still possible, but 
the expressions for the potential in those regions not covered by the 
results of §§ 3, 4 are no longer of simple analytical form, and it would be 
necessary to have recourse to numerical evaluation in general. 

The simplest method of procedure is to determine the function f(«, B) 
of eqn. (17) systematically over the plane z = +0 and then to calculate 
the potential from eqn. (18). A statement of the formulae in full would 
be very lengthy and is easily reproduced, so only a brief statement of 
method will be given. We consider first that part of the problem for 
which ¢, = 0 on z = 0 except on S. On referring back to Fig. 1, the 
potential in the region M'MT7, 7 (for a wing of much smaller aspect 
ratio than is shown) is 


r B Ax(a’) 
‘ # da * p(x’, B’) dp’ yaa .B’) dp’ . [ ) dp" _9 


vie—-B) * J VB-B) 


. (a x) | ‘ (Pp B’) 
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where we have v(a, 8) = 0 when a < ap. Similarly for the region L’' LT, 7’ 

we have : 
B BiB’) Bi(p’) 

* dp’ A l j . 

‘ | dp’ ( » x’, B’) "da! (a’, iets. | p(x B’) da’| 


V(B—B’) Node i oe a V(a—o’) V(a—a’) | 





(54) 
and p(a,B) = 0 when 8 < B,. 


Equations (53) and (54) are satisfied if the contents of the brackets 
vanish (ef. eqn. (27)), and by solving the two resulting equations, we obtain 


(a, B) meets —_ ? 5 ; x 


—. 


A;() 
& (x, B’)W{Ag(a)—B'} dB" [ v(a, B’W{As(a)—B'} dp’] 


: — ——|, (099) 
; 3—p" B—Pp’ 
A;(a) 0 
] 
a 
via p) ami{a— B,(B)} 
B2(B) ; ; B,(B) ; 
A(x’, B),/ ss nd Sasa ea | Ha", Bt Bal (B)— aif d (56) 
= xX a’ - o— x’ 
Bi(B) 0 


When 0 < a < ap, v = 0 and eqn. (55) gives p in this strip (ef. eqn. (28)); 
similarly v is given in 0 < B < B, from eqn. (56). A second application 
of eqn. (55) then gives p in the strip 1) < « < B,(B) since v is now known 
in this strip and so on. Thus p» and v may be calculated systematically 
in the regions for which these equations are valid. 

An alternative procedure is to substitute the right-hand side of eqn. (56) 
in the right-hand side of eqn. (55) and vice versa, thus obtaining integral 
equations for » and vy, which may be solved by the process of successive 
substitution. Since » and v vanish for « < 0 and B < 0 respectively, the 
two series of repeated integrals obtained will terminate in general. The 
series do not terminate, however, if Z and M coincide, that is, if the wing 
is pointed and its leading edges lie inside the Mach cone from its vertex, 
so that the method of this paper is not ideally suited for a treatment 0! 
this case. A solution of eqns. (55) and (56) in closed form would be 
valuable in this connexion, but the author has not been able to achieve 
this. 

In the regions 7, 7, NN’ and 137, PP’, account must be taken of the 
effect of the vortex sheet downstream from the trailing edge. The method 
for doing this is exactly the same in principle as that used already in 
§§ 3, 4, and y, v, p,, v, may be calculated systematically in the appropriate 
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regions in the way indicated above by solving a succession of Abel type 
integral equations. 

The potential may then be calculated everywhere on the wing from 
eqn. (18) and the velocities obtained by differentiation. It must be stated 
that the above process is likely to be extremely tedious in general. 

Once the potential is known all over the wing, the potential may be 
found everywhere else by a process already given by the author (ref. 5). 
The potential in z2 > 0 is given by eqn. (4), the area of integration being 
those portions of S and T for which the inequality (5) holds, since ¢, = 0 
over R. By integrating eqn. (4) by parts with respect to x’, and remember- 


ing that ¢, is continuous on z = +0 and vanishes on the leading edge, 
and that 0d, Ox 0 on T, we have 
we fC Od, (2, y’, | 0) 
0, (%, Y, 2) — Xx 
7 CA 


2(a—x’) dx'dy’ 7” 

(wv +e? By—yp Bey 

the integration now being taken over that part of S alone for which the 

inequality (5) is satisfied. In its general form, this result (57) ranks in 

importance with eqns. (3) and (4), and is often more convenient than 

eqn. (4) in applications, since it is no longer necessary to consider ‘finite 
parts’. 


6. The aerodynamic forces on the wing 

The aerodynamic forces may be calculated by considering the rate of 
change of momentum through a closed surface enclosing the wing. This 
surface may be chosen to consist of a cut tubular surface of circular 
cross-section having a small radius e with the centres of cross-section on C, 
and two plane surfaces lying very close to the upper and lower sides of S 
oining the cut edges of the tube. If n is the outward normal to this 
surface at any point, then the vector force F on the wing is given to our 
order of approximation by 


F ( (pn pU*|(1— M? é¢/ex)i.n-+-n.V¢d|V¢} dS, (58) 


where p is the density of the undisturbed stream, p is the disturbance 
pressure given by the approximate Bernoulli equation, 


} d 22 /D4\2 at\ 2 ad\ 2 
p = —pus(e_ = (PY 4 2 (2) 42 (4) |, (59) 
| Ga 2 \ox 2\ coy 2\ez) } 


and the integration is taken over the enclosing surface. 
If dy) and ¢, are respectively the symmetrical and antisymmetrical 
potentials of § 2, then the contribution from the plane surfaces to the 
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drag force measured in the direction of the undisturbed stream is 
2pU2 v | | ee eects 4 a fy dxdy (60) 
02 02 Janse 
on taking account of the symmetry siiiiai of d, and 4,, the integration 
being taken over the upper plane surface. 

In order to determine the contribution from the tubular surface, we 
introduce locally cylindrical coordinates r,4,s at C: s is distance measured 
along C, r is distance from C, and @ is measured from the upper side of §, 
being zero there and taking the value 27 on the lower side. We require 
the value of the integral over the tubular surfaces as «-> 0. Generally 
Vy) is finite everywhere and |V4,| is finite except on the subsonic leading 
edges, where, from eqns. (55) and (56), it is seen that 04,/ez2 ~ Kr- on 
6 = wand é¢,/éz ~ (constant) on 6 = 0 or 27 as r + 0, K being a function 
of s. If & is the inclination of the leading edge to the x-axis, then the 
contribution from the tubular surface is 


lim pl y2 | ds 
e—-0 a 
T2L,MT 


I< .. M? sin? 2u cos 0 ody Leos A~ ce sin IV's 


or ré | 


where V’ is the vector gradient operator with respect to 2 and y, and the 
> < 


27 


€ d@, (61) 











rPe~€ 


integration with respect to s is taken along the subsonic leading edges. 
All other terms in eqn. (58) vanish in the limit. 
The quantity 
Qn 
od, or lod Ch, . , ” 
lim pU? | ie — M? sins cos 0{— cos 6 — — psiné) |v dy, | e dé (62) 
. 


e—0 cr r< € 


is the dies force per unit length on the subsonic leading edge analogous 
to that occurring in subsonic theory, and is directed normally outwards 
from the wing edge. It is shown in the Appendix that the magnitude of 
this force per unit length is 
apU?K? 
J (1— M2? sin*p) 


The function K may be obtained from the expressions found for p and: 


U2F . (63 


When the flows at the subsonic leading edges are independent, we maj 
use the result for » given by eqn. (28) at the edge M7;, and the similat 
expression for v at the edge LT,. Near the edge M7), we have 
A>(x) 
==0 ~ iB a(x) § a \/{A,(«)- B} 


A;(a) 
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If by is the Mach angle (cot py, B) we have near C in z 0 
B—A,(a) 2r cos ys, cosec(s,,—w), (65) 
A2(x) ‘ ’ 
1B ’ — X(«, B’) dB 


none . {2 cosy, cosec(ypy,—y)} , V{A4(a) p (96) 
Ai(a) 
sy adding the contributions from the plane and tubular surfaces in the 
limit as « > 0, we obtain an expression for the drag force on the wing in 
the form 
Drag 
Lol? 2 OL OZ Ox 


. 


— es 
4 - Po 1 “Pa ) dady—2 F dy, (67) 
i ore T2L,MT; 
the integrations being taken over the upper side of S, and along the sub- 
sonic leading edges respectively. 
A similar expression may be obtained for the side force on the wing in 
the direction of the y-axis. 
The lifting force in the direction of the z-axis is given to our order of 


ypproximation from eqns. (58) and (59) by 


6 lod 
2 ol? | q ) dxdy. (68) 
7 « \ O% c=+0 
Since 6 = 0 on the leading edge, integration with respect to x gives finally 
Lift { 
TIL 4 (d,). o dy, (69) 
oa . 
T2PNT, 


the integration being taken along the trailing edge of S, or, as a result of 
condition (12, ii), along any line spanning T. 
The pitching, rolling, and yawing moments of the aerodynamic forces 


may be found in a similar manner. 


APPENDIX 


In order to evaluate the integral of eqn. (62) for the suction force per unit length 
t the subsonic leading edge, we transform eqn. (1) to a system of locally orthogonal 
oordinates s,t,z at C, where s is distance measured along C, and ¢ is distance from 
C such that ¢ rcos@. In a neighbourhood of C, the variation of ¢, with s is small 
mpared with the variation with ¢ and z, and the equation for 4, is 
od, 9) Y > 2 9 
a4 0, 4 1— M2 sin2us (A. 1) 
ct? O22 
approximately 
Let Cz, t and define new variables 7,, 6, by 
y, 0086, = t, r, sing, = z. (A. 2) 
= — 
Then a solution of eqn. (A. 1) behaving in the required manner (@¢,/¢z ~ Kr-! when 
g Tas? » 0 


2Krt cos 30,-+ O(r). (A. 3) 
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From egns. (A. 2) and (A. 3) 


oy on K cos 36, "4 O(1), (A, 4) 
or Crt lo , 
c K 
fr in oat °° 26, +0(1), (A. 5) 
06 es sec*6,_ 4 
06, C4+tan?6,’ (A. 6) 
whence the integral of eqn. (62) is 
72 [| Od, te ody, od; | 
lim pU? | — M2? sin%ys cos 0°?!) 10 
raat ‘ \ cr ™ a ct ot ben” 
0 
— » UK? [cost 46, e004, (1— (Or 1)eon6,) dh 
: Cc? 74 +-tan?9, 
0 
apU?K? (A.7 


V1 — M2 sin2y)” 
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EDDY DIFFUSION AND EVAPORATION IN FLOW OVER 
AERODYNAMICALLY SMOOTH AND ROUGH SURFACES: 
4 TREATMENT BASED ON LABORATORY LAWS OF 
TURBULENT FLOW WITH SPECIAL REFERENCE TO 
CONDITIONS IN THE LOWER ATMOSPHERE 
By K. L. CALDER 
(Chemical Defence Experimental Establishment, Porton, Wiltshire) 


[Received 5 April 1948] 


SUMMARY 


A theoretical treatment of the problems of eddy diffusion and evaporation is 
leveloped, based on the use of well-established laboratory data for the laws of 
turbulent flow. Expressions which are free from adjustable constants and involve 
only directly measurable quantities are derived for the two-dimensional diffusion 
of matter and evaporation of liquid in flow over aerodynamically smooth and rough 
surfaces. For aerodynamically smooth surfaces the evaporation formulae obtained 
are identical with those found previously by O. G. Sutton using a different treat- 
ment. The latter results, which have been shown to be in good agreement with 
laboratory experiments, may be derived independently as a special case by a more 
direct method involving fewer assumptions. The formulae obtained for aero- 
lynamically rough surfaces are in excellent agreement with observational data on 
the diffusion of smoke in the lower atmosphere under adiabatic conditions, and on 
he evaporation of liquid from the ground surface. 

The treatment may also be applied to the transfer of heat from aerodynamically 
smooth and rough surfaces provided the temperature differences are not too large. 


1, Introduction 

THE present paper makes use of a coefficient of turbulent diffusion derived 
from experimentally established laws for the variation of mean velocity 
with height in fluid flow over aerodynamically smooth and rough surfaces. 
Solutions are provided for the problems of the two-dimensional diffusion 
of matter and evaporation of liquid in turbulent flow over such surfaces, 
in the absence of buoyancy effects. The work is essentially an extension of 
the Prandtl (12, 13) approach to the turbulence problem, although it does 
not involve reference to the latter’s mixing-length theory. In contrast 
to an earlier treatment of similar problems by O. G. Sutton (18a and b), 
no use is made of empirical formulations of the Taylor correlation coeffi- 
cient (21 b), and it is shown that Sutton’s theoretical expressions for 
evaporation from smooth liquid surfaces of finite extent, which have been 
substantiated experimentally, can be deduced as a special case by the 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 2 (1949)] 
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K. L. CALDER 
present method. Provided the temperature differences are not too large 
the present treatment may also be applied to the transfer of heat from 


aerodynamically smooth and rough surfaces. 


2. Aerodynamic considerations 

The existence of two fundamentally different limiting types of turbulent 
boundary layer flow is now well known from experiments with pipes and 
flat plates in the laboratory (e.g. Goldstein, 6). In so-called ‘aero. 
dynamically smooth flow’ the drag exerted on the fluid by the boundary 
surface is a function of the Reynolds number only and, except for a region 
very close to the surface, the distribution of mean velocity is represented 
to a high degree of accuracy by the logarithmic law 


where u = the velocity at distance z from the pipe wall or plate, 
‘~ = 4/(to/p), the so-called ‘shearing stress velocity’, 
T) = shearing stress per unit area of the plate or pipe wall, at the 
transverse section considered, 
p = fluid density, 
v = kinematic viscosity of fluid, 
k = 0-4, the so-called von Karman’s constant. 
On the other hand, for ‘aerodynamically fully rough flow’ the drag is 
independent of the Reynolds number and is proportional to the square 
of the velocity, while the mean velocity distribution is given by 
je = - log z fl {2 
di v ai i) 
where k is again von Karman’s constant and z, (called the roughness 
parameter) depends on the size, shape, orientation, and spacing of the 
roughness elements on the boundary surface. It is evident on physical 
grounds that equation (2) should involve a ‘zero-point displacement 
distance, to provide a datum level above which active turbulent exchange 
first commences and the logarithmic law is valid. Thus modified it becomes 


u | z—d sia 
—_ = - log, ——|}, \o 
Ve = 


the zero-point displacement d being of similar magnitude to the depth 
of the relatively stationary air trapped among the roughness elements of 
the surface. This formula only applies, of course, when z is greater than 
(29+). 


Whether the flow over a given surface is of the aerodynamically smooth 
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or fully rough type depends on the Reynolds number, and smooth flow 
eventually occurs in all cases at sufficiently low values of the latter. 
Laboratory experiments by Schlichting (16) have shown that the criteria 
for fully rough, transitional, or aerodynamically smooth flow can be 


expressed in the form 


, Vy Z 
fully rough flow ——— > 2 
v 
i — a 4 
transitional flow 2-5 > +> 0-13 p , (4) 
Vv | 
;, Uy 2 
aerodynamically smooth flow 0-13 > -*“° 
Vv | 


the roughness parameter z, characteristic of the surface being deduced 
from the velocity profile (3) observed at sufficiently high Reynolds numbers 
with fully rough flow. 

Jecause of the intractability of certain of the subsequent mathematical 
developments of the present paper associated with the use of the loga- 
rithmic velocity distributions (1) and (3), it is necessary to replace these 
laws by approximate power laws. For smooth flow one well-known formula 
of this kind is the Blasius formula (Goldstein, 6) 


u oo 1 
8°74(v,y 2 v)?, 


(5) 
which, except near the boundary surface, is found to agree closely with the 
more accurate distribution law (1), up to v,z/v = 700. For higher values 
of v,2/v, it is found, however, that the }th power law no longer gives 


agreement, and that ith, gth, 


£ 
i 


ith, ete., powers have to be taken, giving 


ormulae of the type 

U/ Vs Q(Ux Zz v)*, (6) 
where g and « depend on the range of v,z/v within which (6) is to agree 
approximately with (1). The approximate power law (6) is used in the 
discussion of aerodynamically smooth flow in the present paper. 

A suitable approximation formula for fully rough flow has previously 
been used by Prandtl and Tollmien (12). Comparison with the power-law 
approximation for smooth flow shows that the appropriate approximation 
lor the velocity distribution (3) of fully rough flow is 


) Sis x 
t a(* “) , (7) 
Ux zo 


where the values of qg’ and « will depend on the range of (z—d)/z, within 
which (7) is to agree with the more accurate formula (3). 
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3. Wind-velocity distribution in the lower atmosphere 

The preceding discussion is of obvious importance when the flow near 
the earth’s surface is considered. Atmospheric flow under adiabatic condi- 
tions may be expected to be comparable with that found in pipe and flat 
plate flow in the laboratory. 

Whether the earth’s surface is aerodynamically smooth or rough can 
be answered in any particular case with a fair degree of certainty. If the 
surface is smooth, then equation (1) may be expected to be valid. The 
‘friction velocity’ v, in this equation may be evaluated from the velocity 
observed at a particular height, and the velocities at other heights may 
then be computed from the equation and compared with the observed 
values. This test has been applied by Rossby (15 b) to certain wind- 
velocity observations over the sea, and he concluded that in moderately 
light winds the sea surface may be aerodynamically smooth.+ The test 
has also been applied by P. A. Sheppard (unpublished) to wind observa- 
tions obtained by Best (1) over very closely cropped grass (a cricket pitch). 
For wind velocities between 1-5 and 4-0 m.sec.—! and adiabatic conditions, 
Sheppard found no agreement between equation (1) and the observed 
wind-velocity distribution. His results are reproduced below in Table 1, 
the velocities being in arbitrary units. 


TABLE | 





z (height above surface,| 25 | 5:0 10 25 |} 50 | 100 200 | 506 
cm.) | 

u (observed) : - | 364 49°0 6370 79°2 90°3 100 I1I‘7 122'5 

u (calculated) , ‘ 63°4 70°3 q7°1 86°1 93°0 | (100) ari Be 








Observed wind-velocity distribution under adiabatic conditions over closely cropped 
grass, compared with theoretical distribution for aerodynamically smooth flow. 


Over the height range investigated, the observed velocity gradient is 
considerably greater than that given by equation (1). The inference is, 
therefore, that a closely cropped grass surface is rough, even at the relatively 
low wind velocities considered, and for greater velocities will, of course, 
remain rough. Any grass surface to be found in practice will therefore 
be aerodynamically rough according to this test. However, the test does 
not suffice to show that the flow in the lower atmosphere is normally of 
the fully rough type, as governed by equation (3), and not of a transitional 
type. To establish completely the validity of equation (3) as applied to 
flow in the lower atmosphere it must be verified that the variation of wind 
with height is of the given logarithmic form and also that the implied 


+ This conclusion has since been confirmed by Montgomery (9) and Sverdrup (20a), 
from observations of the vertical humidity distribution above the ocean surface. 
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relation between surface drag and the wind velocity observed at any 
height, namely 

; . k?u* 
alan (8) 
f ~ ~ 12 
(log.(z—d) /295 
holds in the atmosphere. 

Extensive measurements recently made by E. L. Deacon (unpublished) 
over various grass-covered surfaces on Salisbury Plain have shown con- 
dusively that the wind-velocity profiles in the lowest decametre of the 
atmosphere under adiabatic conditions can be represented with great 
accuracy by a logarithmic law of the form given by equation (3). A brief 
summary of his results (for adiabatic conditions) is given in Table 2. 


TABLE 2 








luerage ratio Zero-point Roughness 
p Wind ve ty of velocities displacement parameter 
(cm.) 1m. (M1.S¢ tom, |/Uym d (cm.) Z (cm.) 
o-7 I 1°45 +15 15°9 
7 1°35 16 8-8 
7 3 1°32 +21 ! 5°6 
7 "5 1°28 + 32 370 
“so I*112 © approx. | 0°20 
1-8 I*I40 ° o-71 
I‘1g! re) ” 2°65 
1-170 ° I°7 





Parameters in the logarithmic law (3), related to grass length, from observations by 
£ £ A 


E. L. Deacon on Salisbury Plain. 


It is seen that over short-grass surfaces the zero-point displacement d 
can effectively be taken as zero. Over long-grass surfaces, however, d is 
definitely positive and becomes more nearly equal to the height of the 
vegetation as the density of the latter increases. Deacon’s observations 
support this concept that d is an ‘aerodynamic length’ which is related 
toa linear dimension of the roughness elements and to the density of their 
packing, 

As to the roughness parameter 2), Best’s (1) short-grass observations 
show no significant change in the velocity gradient, as indicated by the 


tatiO Uy m,/Uym.. With changes of wind velocity, implying a constant value 


of z), independent of wind velocity. The same is true of Deacon’s data 
for short grass 2-3 em. long. Deacon’s values for grass 4-4-5 em. long, 
however, show a tendency for the wind-velocity ratio ug ,/U;m, to decrease 
with increasing wind velocity, implying a decrease in 2, while for long 
grass this effect is accentuated. The explanation of these results is that 
the deflexion of the long grass blades by the wind results in a ‘smoothing’ 
of the surface. Clearly in the case of natural surfaces, with complex 
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roughness elements which may themselves be distorted by the wind, the appro 
only practicable method of estimating the values of the parameters d and | The w 
%) in equation (3) is by observation of the wind velocity at three heights ] the di 
at least. surface 

The verification of equation (8), which is easily done for flow through | toam 
rough tubes in the laboratory, where 7, follows directly from the pressure | K(z) a 
gradient along the axis of the tube, is not nearly so simple in the atmo- 
spheric case, since direct measurements of 7, are not available.t Never- 


theless, Rossby and Montgomery (15a) have shown that this equation, | ,pere 
with a value of z, deduced from observations of the wind profile over the m 
open grassland, gives values for the surface drag 79 in close agreement | jhe y 
with those deduced by Taylor (21a) from an analysis of data of pilot | gom_ 
balloon ascents over Salisbury Plain. Cor 
Knowing the roughness parameter z, of any natural surface as obtained even 
from wind-profile observations using equation (3), it is possible by means rough 
of the Schlichting criteria (4) to obtain information as to the possibility | ogay 
of flow of the smooth or transitional types. repre 
On calculating vy 2)/v from values of vy and z, obtained using equation smoo 
(3), for Deacon’s short- and long-grass data of Table 2, the following values be aj 
are obtained, for an assumed wind velocity of 5 m.sec.-! at 2 m. height. | py 4 
Short grass (1-3 em.'high): d—=—0Ocm., 2 = 0:5 em., vy2/v = 110 . 
Long grass (60-70 em. high): d = 30 em., z = 3-0 em., vy 2/v = 990 Simi 
Since the ‘friction velocity’ v, is directly proportional to the wind velocity, heigl 


this shows that only for the short grass, at wind velocities (at 2 m.) less 
than 0-1 m.sec.-1, can any departure from the fully rough régime be 
expected. Since from measurements of wind profile by Sverdrup (20a) Equ 
the roughness parameter z) even for a smooth snow surface is approxi- 


mately the same as for a short-grass surface, it may be concluded from the wie 

above that, under adiabatic conditions, flow in the lower atmosphere over 
the land is always of the fully rough type. 

, H Sariter 7 an 

4. The eddy diffusivity over smooth and rough surfaces = 
In the turbulent layer over a smooth or rough surface there must be 

a region not too far removed from the surface in which the turbulent Thi 

shearing stress 7, in the mean direction of flow, can be regarded as effec- bot! 
tively independent of distance from the surface. Ertel (5) has shown that 
in the atmosphere the turbulent shearing stress is, to a high degree of 

whe 


| P. A. Sheppard (17a) has measured directly the drag for wind blowing over a very 
smooth extensive concrete surface, but the flow for this rather artificial case was probably T 
of the aerodynamically smooth or transitional type. and 
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ypproximation independent of the height, in the range 1 <z < 30 m.T 
The wind direction is also constant within this layer and coincides with 


the direction of the frictional drag exerted by the wind on the ground 


surface. Above the layer, wind direction changes with height according 
toa modified Ekman spiral. Now the vertical component of eddy viscosity 
K(z) at any height is given (Brunt, 2) by the relation 

T pk ; : (9) 

Cz 

where p is the fluid density and @u/éz the gradient of velocity normal to 
the mean flow. It follows that in the region of constant shearing stress 
the vertical distribution of eddy viscosity can be deduced immediately 
from the vertical distribution of mean velocity. 

Consider first the power-law representations of the velocity profiles 
siven by equation (6) for smooth surface flow and by equation (7) for fully 
rough flow. With the proper choice of constants to fit the more accurate 
logarithmic distributions given by (1) and (3) respectively, these equations 
represent the velocity distributions in the turbulent boundary layers over 
smooth and rough surfaces to a high degree of accuracy. They will now 
be applied in the region of constant shearing stress. From equation (6), 
by rearrangement 

To p(u/q)tA+%(p/z)2u+o9, (10) 
Similarly, by taking a new origin at distance d above the datum plane of 


height measurement for the rough surface, equation (7) gives 


To p i (2) F (11) 
” hee 


Equations (10) and (11) are special cases of a more general form 


7) = epu?B(§/z)208, (12) 
vhere for smooth flow, ) 
3 mS : o = »#, € (1/q)#a+% | 

l+oa ° (13) 
and for fully rough flow, | 
B ® ) Zo © 1/q’?. J 


Chis fact enables a uniform mathematical treatment to be adopted for 
both smooth and rough surfaces. Since in both cases 
u x 
== (z/z,)%, (14) 
1 
where wu, is the velocity at a standard height z,, it is readily verified by 
} For other discussions of this subject see Rossby and Montgomery (15 a), Calder (3), 
and Sheppard (17 b). All assume straight, parallel, steady streamlines. 
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differentiation and substitution for @u/éz in equation (9), using equation 
(12), that 
K(z) = MuP-1z201-2B)z1-a, 


99 15 
where M = (e/x)82%8, U9) 


assuming that in the region of constant shearing stress the shearing stress 
has the surface value 7». 

The above procedure may also be applied to the more accurate loga- 
rithmic velocity distribution laws as given by equations (1) and (3), and 
leads to a simple linear variation of K with height, 

K Vy kz (16) 
for both smooth and rough flow. The use of logarithmic velocity distribu- 
tions, however, leads to insuperable mathematical difficulties in the subse- 
quent treatment of the diffusion problem, as it gives rise to a differential 
equation (see paragraph 5) containing a transcendental function (wu) and 
an algebraic one (K). This type of equation has not yet yielded to treat- 
ment, and consequently the following discussion will be based entirely on 
the use of power laws for the velocity distributions. 

In considering the diffusion of matter, it will be assumed that the 
diffusion coefficient for turbulent mass transfer (eddy diffusivity) is equal 
to that for turbulent momentum transfer (eddy viscosity). Thus if y is 
the concentration of matter per unit volume (measured in g.cm.-*), the 
vertical flux of matter per unit horizontal area, and per unit time, is given 
byt ; 
F= —K—~, (17) 

02 
the negative sign indicating that the flux is directed upwards, i.e. in the 
positive z-direction, when y decreases with height. 

The question now arises as to the validity of the formula (15) for eddy 
viscosity (or diffusivity) for small values of z. Actually over both smooth 
and rough surfaces a transitional region must occur close to the surface, 
in which molecular viscosity and turbulence are of comparable importance. 
In this region, equation (9) should be replaced by 

7 = p(K- vy), (18) 

oz 
where pK(éu/éz) represents the so-called ‘Reynolds apparent shearing 
stress’ and pyvéu/éz is the molecular shearing stress, the latter being 


, - , : , ‘ -o(x , 
+ Strictly speaking the vertical flux is given by F pK —( ), where x/p is the mass 
0z\ p) 


of diffusing substance per unit mass of air. In what follows it will be assumed that the alr 
density is independent of height, so that the flux reduces to the form given in equation (1 i). 
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negligibly small compared with the first except near the boundary. Like- 
wise, near a boundary, the flux equation (17) should be replaced by the 
more accurate equation 


F (K-4 bX, (19) 


Ce 
where D is the molecular diffusivity. Montgomery (9) has actually applied 
this equation to the problem of evaporation from an infinitely extended 
smooth water surface, assuming that the eddy diffusivity K is operative 
down to the top of the laminar sub-layer, where it vanishes. The applica- 
tion of a similar treatment to rough surface flow is, however, not so 
straightforward, and in any case the use of an equation of the form of 
(19) leads to mathematical difficulties in connexion with the main prob- 
lems of the present paper. Various assumptions have been made by 
different writers concerning the nature of the turbulent processes very 
close to a rough surface. Sverdrup (20a and b), considering the problem 
of evaporation from the oceans, assumes that even over a rough surface 
there is a thin layer of air next to the surface, within which diffusion is 
molecular, and applies the laws for fully developed rough flow down to 
the top of this laminar layer. Montgomery (9), however, following Millar 
(8), postulates that the velocity distribution for aerodynamically smooth 
surfaces also applies for a limited distance above a rough surface. Clearly, 
near a rough surface conditions are likely to be complex, and the validity 
of these various assumptions is by no means evident. In the following, 
it will be assumed that the eddy diffusivity, as given by equation (15), is 
valid down to z = 0. The only justification for this procedure is that it 
leads to a simple mathematical treatment and to results which are in very 


good agreement with observation. 


5. Diffusion of matter from an infinite line source 

The problem of the turbulent diffusion of gas or smoke in the lower 
atmosphere, from a continuously emitting line source of infinite length, 
situated at ground level and perpendicular to the mean surface wind 
direction, is one of considerable intrinsic interest. With axes in which Ox 
is downwind, Oy across-wind, Oz vertically upwards, and the source 
defined by (0 f y ©,0) this problem involves the solution of 
the partial differential equation 


Cy c C 
u—X (xe), (20) 
Or C2 Zz 
x being the concentration of matter in g.cm.~ at the point (a, z), subject 
to the boundary conditions 


W) x>VUaszx->oOO, 


5092.6 
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se -o ‘ , a 
(ii) KS +0 as z->0, the ‘ground’ boundary condition of zero 


vertical flux of matter. 
(iii) x is such that the rate of transfer of matter across unit width of 
oc 
a plane x const., i.e. ( ux dz, is constant and equal to Q, the 
0 
constant rate of emission of the source in g.cm.~!sec.-! 

(iv) x > © at the line of emission (x = 0 = 2). 

This problem was first considered by Roberts (14) for the case of 
constant eddy diffusivity. Subsequently he showed, in an unpublished 
paper, that the solution of equation (20) subject to the given boundary 
conditions, and with wind varying as a power of the height w = Uz", and 
eddy diffusivity varying as a power of the height K = Az”, was 


x 


Vl TH(m+1)(m—n +2)}- Vexp{— U: Tym n *®/A(m— = -n , 2)? x} 
(m— n+ 2) 2m +Din—n+2)}—) im +Dlom—n +2) J {(m-+-1)/( (m- —n+2 


atom +1) n+2 
(21) 
which is valid for (m—n-+-2) > 0, I'(@) being the Gamma-function. 
Applying this solution to the particular variations of wu and K given by 
equations (14) and (15) respectively, we obtaint 


Qexp{—ujt- —B)y20+1 M(2 » | 1)2z20 Bhar (22) 
' ae |} (20+) M4 DiQa+ DT{(a-+1) )/(20 FE Djugere a (x-+1)(20+1)? 5 
eee K = (206+ 28—1)/(2«+1). 


This general formula, with suitable values for the constants, again applies 
to both smooth and rough surfaces. 
For smooth surfaces, from (15) and (13), 
M 1 (1 q)2A+@y2al(1+o # (23) 
X 
and (22) reduces to 
Qexp{— uz™t+%220+4/M (2a 4 tae ey a 
X (2 y+] yi (2a4 1) MM X+1)/(2ax DI (ax aie +1)/(2a0+1) )Judlee + Dz; x/(2a4 Val x+1)/(2a+1) 
(24) 





Defining the height H of the cloud of diffusing matter as the distance 

from the surface z = 0, to the point at which the concentration has fallen 
to one-tenth of its surface value, it readily follows from (24) that 

‘ ¢ — 9a l(1 4-0) ~202(1 4-0) 120-4 OF 

H = {2:3M(2a+1)Puz 2la+mgsaild tog M2a+n), (20) 

+ For rough surfaces it should be noted that z 0 is at a distance d above the datum 


plane of the surface, and the source is assumed to be at this distance, where turbulence 


first becomes operative. 
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Similarly for rough surfaces 


M = 22*/aq", (26) 


and (22) reduces to 


Q exp{—2z?%+1/M(2a-+ 1)*2} 


OEE fase eon Bie: = 27 

X ~~ (94+ ] Neat) Ye+viea DP (a+1)/(2a+1)}a, zp art Diea+D’ (27) 
while the height of the cloud is given by 

H {2-3M (20+ 1 )2g-}K2a+0), (28) 


It will be noted that for rough surfaces the vertical distribution of con- 
centration, and consequently the height of the cloud, are independent of 
the wind velocity, while actual concentrations are inversely proportional 
to the wind velocity at a standard height. For both smooth and rough 
surfaces the concentration varies with distance, inversely as a power less 
than unity of the distance, the actual value of which is determined by 
the vertical profile of wind velocity, i.e. the value of « in equation (14). 

To apply the above formulae it is first necessary to estimate suitable 
values for the parameters g and a for smooth surfaces, and q’, «, and Z, 
for rough surfaces. Considering first smooth surface flow, the approximate 
power law (6) must be so chosen that it closely fits the logarithmic law 
1), over the range of v,,2/v with which one is concerned in relation to the 
region of vertical diffusion. The latter proviso is necessary since, owing 
to the approximate nature of the power-law representation, the index a 
isa funetion of the height interval considered, and this means that when 
using a power-law treatment one must strictly have some a priori know- 
ledge of the extent of vertical diffusion for the distance of travel considered. 
lfobservations of wind velocity are available at two well-separated heights, 
which between them include most of the region of vertical diffusion, then 
a suitable value of « can be obtained immediately, using equation (14). 
The value of vy appearing in equation (6) is estimated using the logarithmic 
law (1), from the velocity observed at one definite height, and hence the 
value of g can be obtained. If actual velocity measurements are not 
available at two heights covering the region of vertical diffusion, then 
values may be extrapolated, using (1), and the same procedure adopted. 

Although no observational data are available relating to the diffusion 
oi gas or smoke over aerodynamically smooth surfaces, a rough estimate 
of the vertical diffusion can be made in any particular case from equation 
25), using the values of g and « appropriate to the Blasius formula (5). 
Thus, to consider an example, we find from (25) that for z, = 200 em., 
, = 500 cm.sec.—1, g = 8-74, a = }, and a distance of travel x = 100m., 
the height H of the cloud is approximately 200 em. A more accurate 
estimate of diffusion can now be made by choosing more appropriate 
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values for g and a. Thus, assuming in equation (1) that w(200 em.) = 509 
em.sec.-!, we readily find by trial and error (or, more conveniently. 
using a graph given by Rossby, 15a) that v, = 16-3 cm.sec.-! A suitable 
value for « in equation (14), applying to the velocity profile over the 
greater part of the layer of diffusion, can be obtained from the ratio, 
say u(200 cm.)/u(10 em.), the velocity u(10 em.) being calculated from 
equation (1). It is found that « 0-095, and equation (6) with the above 
value of v, then gives g = 11-9. Thus, in the example considered, for 
diffusion to a distance of 100 m. from the source over a smooth surface 
the approximate power law which closely represents the velocity distribu- 
tion over the height of the cloud is 

u _ 119 (Us io 

Vx Vv J 


This differs considerably from the Blasius one-seventh power-law formula, 


(29) 


given by (5), as is only to be expected, since for the example considered, 
vy, 2z/v ranges from 1-08 x 10 to 2-17 x 10* in the layer 10—200 em., values 
which are much in excess of the critical value v,z/v ~ 700, previously 
quoted as an upper limit for the Blasius formula. Clearly, for greater 
distances of travel from the source, the vertical region of diffusion will be 
greater and the appropriate values of « will be even less. This question 
will be considered further, for the rough surface case. Using the above 
values for q and « in equation (25) leads to a more accurate estimate, 
for the height of the cloud at 100 m., of H = 350 cm. 

Likewise for rough surface flow, if observations of wind velocity are 
available at two well-separated heights which include the region of vertical 
diffusion, then a suitable value of « can be obtained immediately using 
equation (14). The values of v, and z, in equation (7) are estimated using 
the logarithmic law (3), from the velocities observed at three different 
heights. Substitution of the values of v, and « so obtained in (7) then 
gives q’. If actual velocity measurements are not available over the region 
of vertical diffusion, then values must be extrapolated from the logarithmic 
law, and the same procedure adopted. In Table 3 are given values of 4 
and «a so deduced, by extrapolation of the logarithmic law determined 
from the velocities observed at 200, 300, and 500 cm., over a long-grass 
surface, up to the greater heights of 1,000, 2,000, 5,000, and 10,000 em. 
The velocity at 200 cm. was in all cases used as a lower reference height 
velocity. Although the extension of the logarithmic law to the extreme 
heights is probably not completely justified, the values of q’ and « deduced 
are not likely to be in serious error. Table 4 contains similar values as 
determined from the velocities observed at heights of 200, 300, and 500 
cm. over a short-grass surface. 
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TABLE 3 


ght 2 (cm.) 200 500 1,000 | 2,000 | 5,000 10,000 
Observe Extrapolated 
~ aa — Se 
Velocit (cm.s ) 5 55 624 715 503 917 1,002 
¥ ‘220 | 0205 | o194 o'179 o°I70 
$14 4°42 4°61 4°91 5"10 








Values of g’ and a in equation (7), for a long-grass surface, and appropriate to various 


eight ranges (v's 49-5 em.sec.~!, d 30 cm., 2 3 cm.). 


TABLE 4 








g 500 1,000 2,000 5,000 10,000 
( Extrapolated 
' “a i a 
1 ( 5 ) & 52] 574 632 690 766 82 
| 
x 153 o140 | o*140 oO°132 o"125 
6°01 6-28 6°50 6°32 6°97 
Values of g’ and « in equation (7), for a short-grass surface, and appropriate to various 
ght ranges (7 33°3 em.sec 7: 0 cm., 2 0-5 em.). 


“0 


Observations show that the height of a smoke cloud under adiabatic 
conditions over downland is about 1,000 em. at a distance of 100 m. from 
the source. Thus, in this case, the height range 200—1,000 em. covers the 
major portion of the region of vertical diffusion and the appropriate values 
of g’ and a are those given in Table 3 for z 1,000 cm. For greater 
distances of travel it may be assumed as a first approximation that the 
height of the cloud increases linearly with distance from the source, so 
that the remaining values of q’ and a in the table are roughly appropriate 
to distances of travel of 200, 500, and 1,000 m. It is seen that the values 
of q’ and « differ considerably for long-grass and short-grass surfaces and, 
in both cases, « decreases and q’ increases as the vertical depth of the layer 
considered increases. However, for a given type of surface, these varia- 
tions are not so large as to prohibit the use, in approximate calculations, 
i mean values over considerable height ranges, so that precise a priori 
knowledge of the region of vertical diffusion is unnecessary. 

In Table 5 are given values of the ground-level concentration and the 
ieight of the cloud from an infinite line source, calculated by equations 
27) and (28) respectively, for travel over both long- and short-grass 
surfaces. The corresponding theoretical values for an aerodynamically 
smooth surface are also included. 

In agreement with physical expectation, the predicted diffusion over 
i short-grass surface is considerably less, ie. greater concentration and 
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smaller height of cloud, than over the rougher long-grass surface with 


which is associated a higher degree of turbulence. Bearing in mind, how. 
ever, that the two cases considered are extreme ones, changes in the rate 


of diffusion of matter in 


variations of ground-surface roughness, may be expected to be com- 
paratively small and certainly much less than those found to be produced 


the lower atmosphere, produced by normal 


by changes of atmospheric stability. 
Unfortunately no observations suitable for testing the theory are 


available for surfaces of s 


uch small roughness as the short-grass surface 


TABLE 5 





: 
| Long-grass surface 
| 


Rough surfaces (equations (27) and (28)) 
é / 





Smooth surfaces 
(equations (24) and (25)) 











Distance Short-grass surface | 
from Vm = 49°5 cm.sec.-?; Ve = 33°3 cm.sec.?; 
Sour CE @ = 30 cm. 3 d=0cm.} Vy = 16°3 cm.sec.}; 
(m.) Zp = 3 cm. Sy = 0°5 CM. | v = or15 cm.*sec.}; 
10¢ q 4°42 Xs = 3,540 q’ = 6°28 Xs = 5,300 q = 119 Xe = 14,300 
x= 0205 | H = 1,050 x = 0146 H = 810 x= 0095 | H = 350 
1,000 q’ = 5"10 Xs = 410 
xX = 01170 H = 6,800 - 





Calculated diffusion under adiabatic conditions over various surfaces, from an infinite line 


source of strength 1 g.cm.'sec. 


(Xs 


ground concentration in mg.m™*; H 


1, in a wind velocity u(200 em.) = 500 em.sec.“} 


= height of cloud in cm.) 


considered above, since, for the theory to be applicable, a very large area 


of uniform roughness is re 


quired to ensure that turbulence characteristic 


of this roughness is developed up to considerable heights. A large body of 


data has, however, been ec 


Station, Porton, Wiltshire (Sutton, 18c), relating to the diffusion of 


lected at the Chemical Defence Experimental 


smoke and gas under adiabatic conditions, from infinite line sources, over 


downland (grass 20-70 em. high) on Salisbury Plain. These data cover 


distances of travel between 100 and 1,000 m. To compare these observa- 


tional results with the predictions of the present theory as applied to the 
long-grass case above, it may be noted that, in spite of the actual 


dependence of the index « 


on the depth of the layer of diffusion, i.e. the 


distance of travel, the associated variations in the indices of z and x in 
equation (27) are quite small, and for distances of travel between 100 and 
1,000 m. a close approximation to these indices will be obtained if the 


mean of the values of « for these distances is adopted (a 


(0-187). In 


Table 6 are collected the observed and calculated properties of the cloud. 
The agreement between theory and observation is very striking, and 


shows that in the lower atmosphere under adiabatic conditions there 1s 
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TABLE 6 











Property Calculated Observed 
Height of cloud at 1 m. downwind . 1,050 cm. 1,010 cm. 
Independent of | Independent of 
wind velocity wind velocity 
t ground level at 1 . downwind | 3,540 mg.m.-3 3,560 mg.m.*f 
ncentration with wind velocity . | u} u- 
Variation of concentration with distance (100-1,000 0 86 x 0-8 
aR 
Law of variation of concentration in the vertical . | exp(— 21°37) lexp(—2!*) to exp(—z!**) 





Observed and calculated properties of the cloud from an infinite line source under adiabatic 
mditions, over downland on Salisbury Plain. 


Source strength Q 1 g.cm.“'sec.-!; wind velocity at 2 m. height 500 em.sec.~*) 


+ The value quoted is the weighted mean of a number of results, and differs slightly 
from that quoted recently by Sutton (18c). 


close similarity between the diffusion of momentum and matter, and the 
diffusion of matter can be satisfactorily predicted in terms of the vertical 
profile of wind velocity near the ground. 

Returning now to the equation (24) for aerodynamically smooth flow, 
this is identical in form with that derived by Sutton (18 c)—his equation 
(18)—making use of an empirical formulation of the Taylor (21 b) correla- 
tion coefficient and mixing length concepts. Sutton has suggested that 
one reason why this equation predicts a height of cloud at 100 m. of about 
2m., while that observed in the atmosphere is about 10 m., may be that 
the ‘scale’ of atmospheric turbulence is different from that arising in pipe 
flow. In particular, he suggests that in the atmosphere the value of von 
Karman’s constant k may be considerably greater than the value k = 0-4 
as determined in the laboratory. The present development makes it clear 
that this hypothesis is unnecessary, and that the more probable explana- 
tion for the non-applicability of equation (24) to the atmosphere is, as 
ilternatively suggested by Sutton (18 c), that the earth’s surface is not 
aerodynamically smooth. As shown above, the customary laboratory 
value of k — 0-4 leads to results in excellent agreement with observation 
it the rough nature of the flow is appreciated. Further light on this im- 
portant question comes from a consideration of the problem of evaporation. 


6. Evaporation from smooth and rough surfaces of finite extent 
downwind 
The problem of evaporation from a permanently saturated or free liquid 
surface, of finite extent downwind, level with the earth’s surface and of 
such dimensions that it produces no sensible variations of wind structure, 
Was first considered by Sutton (18b) on the basis of his theory of eddy 
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diffusion. It involves, for the two-dimensional case, the solution of the 
differential equation of diffusion (20), subject to suitable boundary condi- 
tions. To formulate the problem mathematically, Sutton assumed that 
for an infinite crosswind saturated strip, of length x), downwind, and 
situated at ground level (z = 0), the appropriate boundary conditions? were 


(i) lim y(x,z) = x, (a constant) (0 <a < ap) 
2—0 


(ii) lim y(#,z) = 0 (0< 2 <2) |\. (30) 
(iii) lim y(x,z) = 0 (0 <2) 


xr—0 

As will be clear from the later discussion, Sutton’s treatment is restricted 
to aerodynamically smooth surfaces, and for the latter it seems reasonable 
to assume that, at the evaporating surface, the constant concentration 
x; can be identified with the saturation vapour concentration of the liquid. 
In the present paper the additional problem of evaporation from aero- 
dynamically rough surfaces will also be considered, such as evaporation 
from an area of the earth’s surface sprinkled with a volatile liquid, and 
therefore not necessarily saturated. For this case, the appropriate surface 
boundary condition is not so obvious physically, as the liquid surface will 
in general be discontinuous, and the area and geometrical disposition of 
the liquid surface will be determined by the quantity of liquid dispersed, 
the method of dispersion, and the characteristics of the surface. As a 
working hypothesis it will be assumed that the concentration of vapour 
over a rough surface, which is saturated or uniformly sprinkled with 
evaporating liquid, is constant over a plane at a distance from the datum 
plane equal to the zero-point displacement d of the velocity profile (3). 
At this distance turbulence can be regarded as first becoming operative. 
Thus the boundary conditions (30) will be assumed valid for rough surfaces 
also, the plane z = 0 in this case, however, being at distance d from the 
datum plane of height reference in (3), and x, not necessarily being equal 
to the saturation vapour concentration corresponding to the temperature 
of the liquid. 

As for the problem of diffusion from an infinite line source of gas or 

+ Sutton considers the case of dry (or vapour-free) air blowing over an evaporating strip 
If the boundary conditions (30) (ii) and (iii) are replaced respectively by 

(ii) lim x(2,z)= x9 (O< 2% < a), 
200 
(iii) lim yx(x%,z) = x9 (0 < 2), 
xr—0 


where xo is the concentration (assumed uniform) in the air stream upwind of the strip, all 
the subsequently derived equations apply, provided ,(x,z) is replaced everywhere by 
x(x,2)— xo, and xz is replaced by xs— xp. 
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smoke, a uniform mathematical method can be employed for both smooth 
and rough surfaces. Using the general expression (15) for eddy diffusivity, 
the diffusion equation (20) becomes 

u,(z/z,)% oe wns I Mf u28—1z20-28)z1—a OX), (31) 

ox ~=—s Gz | dz 

\ detailed discussion of this type of equation, subject to various boundary 
conditions, has been given by W.G. L. Sutton (19), regarding the equation 
as a generalization of the classical equation of heat conduction, and recently 
by Jaeger (7), using the method of the Laplace transformation. In the 
present paper the method of sources will be used, which leads to a simple 
integral equation. 

Considering the infinite crosswind strip of evaporating surface defined 
byz = 0,0 <2 < ap, this may be regarded as an aggregate of elemental, 
continuous, infinite, crosswind, line sources of vapour. If Q(x) dz g.sec.—! 
is the strength of the line source of width dz at x (0 <a < ay), then Q(x) 
is a function of position to be determined by the boundary condition 
30) (i). The general formula for the concentration of matter produced by 
an infinite line source has already been derived in equation (22), which 
is applicable to both smooth and rough flow, if appropriate values are 
taken for the parameters, from equations (13). For simplicity (22) will 
be written in the form 


¢ 
x exp! -y/2), (32) 
a 


where Oo Re 


and for smooth surfaces, by (24), 





N (20-4 ] )U(2a Dyreraea DP 


s+ 1/(2%4+1)4— af(2a+4+1) 
u rd J 9 
1 1 
r+] 


2 
| 
us 1 Kyo Ox l | is 
M (2a-+-1)?290"4 x)? 3 
— | 
while for rough surfaces, by (27), | 
, 11 | 
9 | 1/(2a+1 Y4+1I)(2a+-)l xX ii ? . 
A (20+ 1) M (Sa ns : | 
e2o 1 | 
M(2« 1)? } 
The boundary condition (30) (i) then gives 
6 
V( I \dx . 
Xs (0 <8 < %), (34) 


N(é—2x)? 
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which is a simple integral equation of the Abel type (Whittaker and 
Watson, 22) determining Q(x). Its solution is 


Q(x) = 








sin ond ? Nx, 40 
(x—6)1-2 


© 


a7 dx 
0 


, soz 


— VX5 axel, (35) 


— 
7 





The total rate of evaporation EF per unit width crosswind of the evaporat- 
ing strip is thus 











E (29) = [ Q(a) dx = Ny, 279g, (36) 
0 
while the vapour concentration at a point over the strip is 
fr Q(9) (—y | 
(x )= | —————-_~ OX, —-) dH 
x N(x—6)° P\ — | 
0 
— 6°-1(4— a)-°exp| —% d@ (0<a<a,) (37) 
0 
and at any point downwind beyond the strip is 
sinom [- t af i 
x(z,2) = x, — 6° alla d0 (x > a). (38) 


Clearly these expressions for y(x,z) also satisfy the boundary conditions 
(30) (ii) and (iii). They have also been obtained by W. G. L. Sutton (19) 
directly by solving equation (31), subject to the given boundary conditions. 

The expression (37) may be further simplified to give an Incomplete 
Gamma-function, since the integral is a degenerate case of the confluent 
hypergeometric function. Thus writing @/(~—6) = v and p = y/z, in the 
integral, J say, of equation (37), then 


d0/(~a—@) = dv/(1+v), 


and I pI—le—p(l+ i) We . 
a 1+ v 
0 
dl r 
Hence —_-=>— | yo—le-pll+r) dy 
dp } 
0 
@o 
oe? | s°-le-§ds (s = pv) 
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= —I(c)p-’e-?. 
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r and 20 
Therefore I (co) t-%¢ t dt. since I — 0 for = 0, i.e. for p= 00, or 
p 
I C(o)fT —a)—I(p, 1- o)|, 
where the second term in the square brackets is the Incomplete Gamma- 
function, defined by 





(35 
g 
0rat- 1'(0, nN) {™-1e-t dt. 
0 
Hence from (37) 
(36 sin om * - " 
x(x, z) Xs P(o)[ T(1—o)—T(y/x, 1—o)| 
sin o7 ,, . ‘ 
x.) : — I(o)T(y/x, 1 —o)| (0O<4 <a). (39) 
While this mathematical reduction is possible for expression (37), the 
integral occurring in (38) for the concentration downwind of the strip has 
(37 not so far been reduced to known functions and must be evaluated by 
numerical integration. 
From the above general expressions for H and y we find, using the 
relations (33), that for 
(38 (i) Smooth surfaces (M given by equation (23)): 
. Zat+lly, .. {(a+1)z7 " . +1] 
. E(x) \Xs cin | < as (20-4 ] M2a+) Y(a+Di2e+DP fsck. Berd x 
tions latlfa \(2a+1)} Qa+1 
(19 1f(2a4+1)+— af(2x+1)>(44+1)(2«+1) 
ions x Uy a4 ea ee (40) 
let l . ((a+1)7)./a+1\, y2ul+a)z2a+1 x 
p x(x, Z) y¥,| 1—- sin! - - I] ~~ rl. : — a 
uent TT | (2a 1)} 2a+1 | M (2a i 1 )2zF0+a)y 2a+1) 
t 
. the 0<x< 2%), (41) 
0 


fs wee l)z) ff . 

x(%, 2) Ke ain) fa x(2a+1)( 9) (a+ D(2a+1) y 
0 

ugva + Ot)ay2 ax +1 


e d@ (x2 >). (42) 


exp! a a 
| M(2a+ 1)22207/0+0( 7 —)J 


(ul) Rough surfaces (M given by equation (26)): 


E(x.) = (= ~) X¢sin{ FY (yt 1) HGH JY (0+ DI2a+4) > 
x ] 7 (2a 1)j 


rs) 27 tafetDiea+, (43) 
2Q0-7 
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1. {(a+1) 77 ] sf z2a+l x | 
Ax,2) = x,| 1—- sin] S27 a) = ea a 
x(@, 2) x, zn NQ orn rn i) \M(2a-+1)2x il | 


(0 <X%< Mp), (44) 
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:, 2) = xe ain (et! 7 { Q-Al2x+1( gr - 9) AX+DI2x+D) »¢ 


\(2a-+ 


{ _z2a+] ) | 
| M(2a+1)%(a—6)f dd (a > Ly). (45) 


For computational purposes equations (41) and (44) are most con- 


< exp 


veniently re-expressed in terms of Karl Pearson’s function /(X,>p),; 
values of which have been tabulated (Pearson, 11). Thus in place of (44) 
we have 

eee —(e+1)| 


: : 46 
Xs | M(2a+-1)2& (2a+1) J () 


Since the integrands of (42) and (45) are infinite at 6 = 0, the numerical 
integrations must be carried out in two parts. Thus in the subsequent 
calculation of the vapour distribution downwind of a rough evaporating 
strip, the integral in (45) was evaluated for 0 <@ < 2,/10, i.e. with @ 
negligibly small compared with x, so that it becomes approximately 


9/10 
1 | in a 0 ] 
ae t+Di2a+) exp i —a(2a+1) gg 
\M(2«+1)x} J 
0 


and the interval x)/10 < 6 < a, was treated numerically. 

The expressions (41) and (44) involve x and z in the combination 
z20+1/y only, from which it follows that for given a, the contours of 
constant y are parabolas of degree (2a+-1), with the surface as axis, and 
the leading edge of the evaporating strip as vertex. Furthermore, the 
distribution of y with z, when z is small is y— x, o 2%, since when ¢ is 


p 
small, [(¢, p) = ll Hence, from (39), when z is small, 
P 
" l-o 
‘ git | (co) Y (47) 
7 (i oa) \x 
On eliminating x between this equation and equation (35) one obtains for 
the local rate of evaporation from an infinitely extended surface 
. (l—o) 1 
( N(y.—y)— es 
v ax Io) yo 


+ Pearson tabulates I(u,p), where u = 
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which can otherwise be obtained by integrating the simplified form of 
diffusion equation, viz. 

“x ey) <0, 

C2 | Oz 
which is applicable when the distribution of vapour can be regarded as 
horizontally homogeneous, i.e. not varying with distance in the mean 
direction of the wind. 

Equation (41), for the vapour concentration above a smooth saturated 
strip, is identical with that first obtained by O. G. Sutton, in an un- 
published paper, by solving the equation of diffusion in relation to his theory 
of turbulence and subject to the boundary conditions (30). The expression 
(40), for the total rate of evaporation from a smooth evaporating strip, 
was first obtained by the present writer, by direct integration from 


equation (41), using the relations 
: { rR , : C 
E(x») (Wx) p=7_ 42 = | (tim Ke) dx. (48) 


It has been shown by Pasquill (10) to be in very good agreement with the 
results of wind-tunnel experiments involving a variety of liquids. - The 
present formulation shows, however, that the final expressions may be 
derived independently as a special case without making use of Sutton’s 
empirical formulation of the Taylor correlation coefficient or of the semi- 
empirical concepts of the Prandtl mixing-length theory of turbulence. 

With a one-seventh power-law velocity distribution, which is valid at 
low Reynolds numbers (equation (5)), equation (40) gives 

Biz.) oc a™aa™, 
an experimental result which has been quoted previously for smooth 
surfaces (Brunt, 2) and which also represents very well the heat transfer 
by forced convection from a heated smooth plate (Elias, 4). 

Turning now to the case of rough surfaces, equations (44) and (45) show 
that, in this case, the distribution of vapour concentration is independent 
if the wind velocity and depends only on the roughness of the surface, 
which determines the degree of turbulence present. Furthermore, equation 
43) gives 

Bi Lo) OC Uy T+ Ue +1) 
80 that for rough surfaces the total rate of evaporation should be directly 
proportional to the wind velocity u, at the reference height z,. Since the 
index of x, is of the same form as that for smooth surfaces, and only changes 
slowly with a, the variation of E with 2, will be approximately the same 


us for smooth surfaces 
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As far as is known, no laboratory observations at present exist for rates 
of evaporation from aerodynamically rough surfaces. Field experiments 
have, however, been carried out at the Chemical Defence Experimental 
Station, Porton, Wiltshire, to obtain information on the distribution of 
vapour concentration both over and downwind of rectangular grassland 
areas (grass length less than 1 inch) sprinkled with aniline under adiabatic 
conditions, the vapour concentrations being determined chemically. The 
contaminated areas were of sufficient extent across wind for them to be 
regarded as of infinite width with respect to the sampling positions used 
on the centre lines of the areas, subject to prescribed variations in wind 
direction. 

Comparison between theory and experiment cannot be made directly, 
since the theory gives the concentration y(2,z) at any point in terms of 
Xs the ground concentration on the contaminated area, and this is not 
known. Satisfactory comparison can, however, be made by expressing all 
concentrations, both theoretical and observed, in terms of the concentra- 
tion at a particular position. 

Tables 7 and 8 give some mean results for concentrations over and 
downwind, respectively, of contaminated downland areas. The values of 
a and q’ selected in each case are appropriate to the average height of 
the vapour cloud. 


TABLE 7 





Distance from upwind edge 











10 yds. 20 yds. 30 yds. 
Height 1” 12: (10°0) r” 11-7 (11°0) r” 12:8 (11-8) 
6”. 5:9 = (47) 9” 5:4 (4°9) | YY 5°5 (4°90) 
” 8 3°2 «6(3°9) 18” 30 (2:2 | 2’ «86390 = (24) 
2 1:0 (1°0) 3 «= to (T°) 4" 0 (10) 





<a canines 

Relative vapour concentrations under adiabatic conditions over an area of short grass- 
land, 30 yds. long, in the mean wind direction and initially contaminated with aniline to 
a density of 50 g.m.~? 


(Calculated values from equation (44): observed values in brackets. 


d = 0-0 cm., z 0-5 em., « 0-25, q 3-56.) 


The general agreement between theory and observation can be con- 
sidered as very satisfactory for field experiments of this type. 

From the above analysis of the two-dimensional evaporation problem 
it may be concluded that the present formulation of the laws of turbulent 
transfer satisfactorily predicts both evaporation from aerodynamically 
smooth surfaces in wind-tunnel experiments and also the distribution of 
concentration from very much larger aerodynamically rough surfaces in 
an adiabatic atmosphere. Since equation (44) or (45) may be used to 
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obtain a value for x,, from a knowledge of the value of x(x,z) at some 
definite position, it should also be possible to calculate total rates of 
evaporation from rough surfaces, using equation (43). Wind-tunnel 


experiments are in progress to test this suggestion. 


TABLE 8 





Distance from upwind edge 











"34 (0°31) | 0°45 (0°83) 


10 yd 1s yds. | 20 yds. 25 yds. 40 yds. 
I I (1-00) 1°00 (1:00) 1°00 (1°00) 100 (1:00) | 1°00 (1:00) 
yutior : "50 = (0°31) ‘91 (0°94) | 0°99 (0°99) I'ro (1°08) 108 (1°13) 
0:79 (0°80) | o89 (1°05) 
o 


013 (0°08) 0°43 (0°51) 0°63 (0°69) 
5s (0°06) o15 (o-r1) 





Horizontal : I'oo (I "62 (0°54) 0°46 (0°47) 38 (0°35) 
distributio1 


0°23 (0°19) 


Q 





Relative vapour concentrations under adiabatic conditions, downwind of an area of 
hort grassland, 10 yds. long in the mean wind direction, and initially contaminated with 


niline to a density of 10 g.m. 
(Calculated values from equation (45); observed values in brackets. 
d 0-0 em., 2 0-5 em., a 0-22, q’ 4-1.) 

The writer is indebted to Professors O. G. Sutton and P. A. Sheppard 
for encouragement and discussion throughout the entire course of develop- 
ment of the work described in the present paper. His thanks are also 
due to Mr. E. L. Deacon, of the Meteorological Office, for permission to 
employ observational results which have not yet been published in full. 
Finally, acknowledgement is made to the Chief Scientist, Ministry of 
Supply, for permission to submit this paper for publication. 
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BENDING OF AN ELLIPTIC PLATE WITH A 


CONFOCAL HOLE 
By B. R. SETH (Hindu College, Delhi) 
[Received 31 August 1948] 
SUMMARY 

An approximate solution for the bending of a thin elliptic plate with a confocal 
hole subjected to uniform pressure and clamped at the edges is discussed. The 
numerical results obtained are compared with those for a complete plate. It is 
found that the maximum deflexion w,, occurs near the inner boundary. For the 
plate bounded by ellipses whose semi-axes are (1-6c, 1-249c) and (1-14c, 0-548c) it 
s found that w,, is almost one-fifth the value for the plate complete up to the 
ter boundary. By making the minor axis of the hole vanishingly small the 
interesting case of an elliptic plate with a clamped crack can be treated. 
Introduction 

TuE effect of a hole or a group of holes on the distribution of stress in 
any member or structure under load is of great importance in engineering 
problems. Of those problems involving forces, isolated or distributed, 
normal to the plane of the plate, there are a number of well-known 
wlutions for plates of various shapes.t Among recent papers extending 
our knowledge of particular solutions of these problems are those of 
Sen (1), Stevenson (2), and Sengupta (3). The finite circular plate with 
an eccentric hole and bent by a concentrated load at any point has been 
liscussed by Kudriavtzev (4). Very recently Southwell (5) has used his 
nethod of relaxation to get numerical results for a plate with perforations. 
It will therefore be not without interest to discuss the bending of an 
elliptic plate with a confocal hole bent by uniform pressure and clamped 
itthe edges. For the complete elliptic plate the solution is due to Bryan (6). 
\s he did not solve this problem directly by using elliptic coordinates (and 
is we shall require some of those results for the sake of comparison with 
those for a plate with a confocal hole), we first consider the complete plate 
problem. By making the minor axis of the hole vanishingly small we can 
liseuss the case of an elliptic plate with a clamped crack coincident with 


the straight line joining the foci. 


Complete elliptic plate 





| The deflexion w of the middle surface satisfies the differential equation 
2 2\2 : 
Viw | le w I . (1 
L" da* * oy? D 


being the uniform pressure and D the flexural rigidity. 


| See, for example, Timoshenko, Plates and Shells, Chapter VII (1940). 


(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 2 (1949)] 
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The elliptic coordinates €, are given by 
a+iy = ccosh(é+7n), (2) 
and a particular solution of (1) is 


l p p ; . 
-L_g2y? — + _ ¢ginh? 2é sin? 2 
sD’ ~ 128D ne 


k| — cosh 4€ cos 4n-+ cosh 4€+- cos 4n— 1], (3) 
where k = pc*/512D. 


The typical solutions of (1) in elliptic coordinates are (7) 


(cos cos 
tenant | 2» (n—1)y+e ne| ~ (n+ 1)y (4.1) 
\sin sin 
cos COs , 
and efen ve ~ (n- 1)y-+4 ceanont| : *(n+1)n. (4.2) 
\sin sin 


In addition we have solutions of the type 
AE+B and ré = hc*€(cosh 2€+-cos 2n). (4.3) 
As w should be even both with respect to € and » we see from (2) that 
w = k| —cosh 4€ cos 4n-+- (cosh 4€—cosh 4x) -+cos 4 |+ 
+-A,|cosh 4€ cos 2n-+-cosh 2€ cos 44] + 
+A,| cosh 2é—cosh 20+ cos 2n|+-A, cosh 2£ cos 2y+ 
+A, cosh 4€ cos 4. (5 
The boundary conditions are 
w=0, dw/o&=—0, for €=a. 


These give rise to five equations involving the constants A;. But it 1s 


found that only four are independent. Thus we get a closed solution fo" 
this case. We find 


4k cosh 2a 4k cosh 2« 
A, = = 3 As -4k cosh 2a; 
1+ 2 cosh? 2a 2-+- cosh 4x ? 


A, ; 16k cosh? 2a, A, 


2k(1-++- cosh? 2x) 
2-+-cosh 4x ’ ; 


2+ cosh 4a 


After some simplification (5) reduces to 


l p zt y\2//3 3 2 - 
) = So | isis: Gian Slade, a ia ——|, (/ } 
” 8 al a ‘| [(Gtita) 
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a.) being the axes of the ellipse. The maximum deflexion w,, is (6) 


lp 3 3 2 
Ww + a on —_ 
= 8 D/ \a* © bt ¢ a®b? 


] p c* cosh*a sinha 





8 D 3(coshta-+sinh‘a)+ 2 cosh2a sinh?’ 


9 
which gives (w,,/k) cosh 4a—4-+ 


m 


2+-cosh 40° 
Elliptic plate with a confocal hole 
In this case we can write with the help of (4) 
k(cosh 4€+-cos 4n)+-A,y &+ Bo 
-[A, sinh 2(¢—8)+ A, sinh 2(—«)—(A, sinh 28+ A, sinh 2a)cos 2y]+ 
| 4, &(cosh 2é-+- cos 2n)+-[A, cosh 2(€—f)-+ A; cosh 2(€—«) |cos 2n + 
| Af cosh(2—a—f)cos 47+ cosh(4é—a—Bf)cos 2n]- 
| A, cosh(4é—2a— 28)cos 4+ 
A, cosh(6£—3a—38)cos 4n-+cosh(4€—3a— 38)cos 6n|-+4 
(A,+ B,é+ A, sinh 2(€—f)+ A, sinh 2(€—a)-4 
+ A,écosh 2€+-k cosh 4€|+-| —(A, sinh 28+ A, sinh 2a)-+ A, &+ 
A, cosh 2(é—f)+ A; cosh 2(é—«)+ A, cosh(4é—a—B) |cos 2-4 
[ A, cosh(2é—a—f)+- A, cosh(4é— 2a— 28)+ 
A, cosh(6£— 3a—38)-+-k|cos 47-4 
(9) 
The boundary conditions are w = 0, éw/é€ = 0 over = a and € = f.. 
[he first term in (9) gives 
1, + Bya+A,sinh 2(a—f)+A, «cosh 2a+k cosh 4a = 0, 
A,+ B,B+A,sinh 2(8B—a«)+A,f cosh 28+-k cosh 48 = 0, 
| B\+2A, cosh 2(a—B)4+2A,+ A,(2asinh 2a+ cosh 2a)+- 4k cosh 4a = 0, 


B,+2A,+2A, cosh 2(8—«)+ A,(28 sinh 28+ cosh 28)+ 4k cosh 48 = 0. 


0 


The second term gives 
1, cosh(3x—B)-++A,+A, cosh 2(a—B)+A3 a—A, sinh 28— A, sinh 2a = 0, 
1, cosh(38—«)+A; cosh 2(8—a«)+A,+A,B—A, sinh 28—A, sinh 2a = 0, 

1A, sinh(3a—f)+2A,sinh 2(a—f)+A, = 0, 


14, sinh(38—.«a)+ 2A, sinh 2(8B—a)+A, = 0. 


(11) 














180 B. R. SETH 





The third term gives 
A, cosh(a—f)+ A,cosh 2(a—8)+ A, cosh 3(a—B)+k = 0, 
A,+4A,cosh(a—f)+3A,[3—4sinh?(a—f)] = 0. (12) 
Solving (10), (11), and (12) we get 
2 sinh(a—f)cosh(«-+-8)—(a—f)[1-+-2 cosh 2(a—B)] 








A,+A, = 2A 
atAs 6 («—B)cosh(a—f)—sinh(a—B) , (13.1) 
A,—A, = 22 sinh(«-+-8)cosh(«—f) 
4 A. Ac —Acosh(a—f)—sinh(a—B)’ (13.2) 
he we hh, sinh?(a—f)sinh(«-+ B) - 





(a—B)cosh(a—f)—sinh(a—f)’ 
2A, sinh 28-+-2A,sinh 2x = A,;—2A,cosh(«+8)[2+ cosh 2(a—B)], (13.4) 
2(A,+A,)[{1+-cosh 2(a—)}(a—B)—sinh 2(a—B)]+ 
+-A,| (a—f)(cosh 2a-+ cosh 28)-+ (a—)(2« sinh 20+ 28 sinh 28)— 
— 2a cosh 2«-+-28 cosh 28]+- 
+-k[4(«—B)(sinh 40-+-sinh 48)— 2(cosh 4a—cosh 48)] = 0, (13.5) 
4(A,—A,)sinh*(a—)-+A,| cosh 21+ 20 sinh 21—cosh 28—28 sinh 28|+ 
+4k(sinh 4a.—sinh 48) = 0. (13.6) 


Combining (13) with (12) we get A,, A,, As, Ay, A;, Ag, Az, and Ay. 
Substituting these values in (10) we find A, and B,. 

It is found from the expression for w given by (9) that it is in general 
quite sufficient to stop at the coefficient of cos4y; if required, further 
terms can be taken. 

To get a numerical idea of the solution we take 


x = 407 = 1-04720, B = 4m = 0-52360, 
and we find 
A, = 2-666k, By = 50-371k, A, = —67-109k, A, = 22-850k, 
A, = —1-074k, A,=0-681k, A; = 0-330k, A, = 0-035k, } (14) 
A,=1:795k, A, = —1-671k. 


The maximum value, w,,, of w occurs on the major axis and its value 
is given approximately by 


4 x 190} 15) 
Wn 5-712k. ( 


The following table gives the values of (w/k) for varying values of é. 





| | | 
E | 05236 | 06 | 065 | o70 | 075 | 080 | ogo | 1°0472 

















w/k o 5°624 | 5°712 | 5°577 | 5:192 | 4°542 | 3°243 
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From (8.1), if the plate has no hole, is bounded by the ellipse € = 0-5236, 
and is clamped there, we have 


(w,/k) = 1-880. (16.1) 
For a clamped elliptic plate bounded by € = 1-0472 
(w,,/k) = 29-259. (16.2) 


For the plate with a confocal hole we see that the maximum deflexion 
yours near the inner boundary. In the present example it is about three 
times the values given by (16.1), and one-fifth the value given by (16.2). 
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PRINCIPLES AND PROGRESS IN THE CONSTRUCTION 
OF HIGH-SPEED DIGITAL COMPUTERS 


By ANDREW D. BOOTH and KATHLEEN H. V. BRITTEN 
(Birkbeck College, University of London, and British Rubber Producers’ 
Research Association, Welwyn Garden City, Herts.) 


[Received 5 November 1947; revised 2 July 1948] 


SUMMARY 

The basic principles underlying the mathematical design of high-speed digital 
computers are discussed and the necessary components of such machines defined. 
Seale of notation, the form of the ‘memory’, the action of the control, and othe: 
practical details are considered, and a brief discussion is included of the exact 
arithmetic functions of which these machines must be capable. 

This is followed by a description of current computer projects in America. Aiken’s 
second relay computer at Harvard, the Bell relay machine, E.D.V.A.C., and th 
Princeton electronic computer are described briefly and an idea given of their stat 
of completion in 1947. 

Introduction 

It is the purpose of this paper to discuss the basic principles underlying 
the mathematical design of high-speed digital computers, and to follow 
this with an account of the progress which has been made, in the U.S.A. 
and at home, with the actual construction of such instruments. 


Digital versus analogue computers 

At the outset it is perhaps worth noticing the essential difference 
between the two possible types of computing machine, digital and analogue. 

An analogue machine is one which performs its functions by finding 
some mechanical, hydraulic, or electrical mechanism whose action in effect 
duplicates the mathematical processes required. A single example suffices; 
suppose that it is required to calculate the function acos@. A simple 
mechanism for doing this is shown in Fig. 1. 

A cross-head OA is capable of rotation about O, a peg B is positioned on 
OA,so that OB = a. BC isa bar rigidly joined perpendicular to CP which 
is constrained to move in the line OCPX. If BC is maintained in contact 
with the peg B and OA is rotated through the required angle @, the dis- 
placement of P from its position when @ = 47 is equal to acos@. It is 
apparent that such a device will have severe limitations in accuracy and, 
for more complicated systems than that considered, an upper limit of 1 in 
1,000 is expensive to realize, whilst 1 in 10,000 is on the extreme border 
of practicability. To paraphrase a remark of Professor Hartree, each 
decimal place multiplies the cost by a factor of about 10. 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 2 (1949)] 
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The most ambitious analogue computer constructed to date (and it now 
seems for all time) is the Bush—Caldwell differential analyser at the Centre 
of Analysis in the Massachusetts Institute of Technology; this achieves, 
on favourable occasions, the upper limit of accuracy quoted above and has 
, tape-feed device which enables a problem to be set into the machine from 


A 
He 





tua tod 


3 
OQ 





Fic. 1. 
a standard commercial electrical typewriter. Electrical gear-boxes elimi- 
nate the laborious setting operation required in previous analysers, and 
the circuits are so arranged that all the integrators receive a fair share of 
the work of the machine. This feature is important in a machine whose 
maximum capacity exceeds that normally required on routine problems. 

The subject of analogue computers will not be discussed further as the 
field is well known and is well covered in existing literature (1). 

A digital computer is defined to be any machine which takes its informa- 
tion in the form of actual numbers (i.e. assemblages of digits) and operates 
with it in this form throughout. It is at once apparent, from this, that only 
the ordinary operations of arithmetic are possible to such a machine, and 
that any continuous process must be replaced by a suitable finite difference 
or other numerical equivalent. The most common example of a digital 
computer is the ordinary desk Marchant, Monroe, or Brunsviga, and it is 
perhaps worth observing that such a machine, with a 6x 6 decimal multi- 
plication time of around 10 seconds, is faster by a factor of about 25 than 
) non-expert (i.e. non-professional computer) using pencil and paper 
methods. At these speeds it is just possible for the human operator to 
supply the machine with data and instructions at an adequate rate. With 
an increase in speed by even a factor of 10 this is no longer true and, if the 
machine is to be hand-controlled, no further improvement in performance 
is justified 


The control of high-speed computers 
These conditions lead naturally to the idea of a machine capable of con- 
trolling its own operations and exerting a certain amount of independent 
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judgement. Before being more explicit it may be well to clear up a point 
which is often raised in this connexion. If a human operator has to 
‘programme’, i.e. instruct the machine at some stage, does this not 
impose a fundamental limitation on the speed? The answer is in the 
negative, since most problems involve a large amount of repetition, and 
if the machine can exercise a little judgement it is sufficient to detail the 
process once and then tell the machine to repeat it, either a fixed number 
of times, or until some criterion is realized. As examples consider three 
typical problems: 

(1) Tabulation of f(#) for « = 1,...,n; 

(2) Iterative evaluation of a}; 

(3) Multiplication of two matrices A. B. 

In the first problem the detailed programme for the calculation of f(x) 
for one value of 2 would be written down. The next value of x would now 
be inserted into the scheme and the process repeated. The operation could 
be terminated by having the machine calculate y = (a—n) at each stage 
and then using a ‘judgement’ order of the type: 

if y > 0 perform operation sequence A; 
but if y < 0 perform operation sequence B. 
In this example sequence B would be that for calculating f(x), whereas A 
would signal the end of the calculation, since (z—n) < 0 until 2 =n. 
The process would therefore terminate automatically after the correct 
number of values of f(z) had been generated. 

The word ‘programme’ as used in this context may require explanation. 
Before a problem is ready for solution in a computing machine it must first 
be broken up into processes which the machine is capable of performing. 
Thus, as mentioned above, when using a digital machine the continuous 
process of differentiation must be replaced by a finite difference formula. 
This translating of a problem in terms of the available functions of the 
machine is generally known as programming. 

The second problem, the evaluation of at to a preassigned accuracy, is 
dealt with in the same fashion although the required number of iterations 
of the process 
Lnay = 34(x,+0/z,,) 


n 


may not be known ab initio. The decision to end the calculation is made 
by evaluating (x,,,,—2,) at each stage and noticing that it is negative up 
to a certain point and zero afterwards (to a preassigned number of decimal 
places). The point at which the two successive approximations become 
equal is the required termination, and a ‘judgement’ order is adequate to 
determine it without human intervention. 
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The third example, matrix multiplication, is effected by a straight- 
forward triple iteration and the same sort of procedure is adequate. 


The scale of notation 


So far in this discussion of digital computers ‘numbers’ and ‘digits’ 
have been mentioned without definition, but one of the first problems 
confronting a designer of machines is the choice of a scale of notation for 
hisnumbers. The decimal, or scale of 10, notation, although hallowed by 
tradition, is by no means mandatory. With present electrical and electronic 
facilities it is an unfortunate fact that most simple devices favour not 
decimal but binary scale, and it is pertinent to consider what is in effect 
the best possible scale. To state the problem in mathematical terms: to 
record numbers up to a magnitude K, what scale of relationship minimizes 
the total number of digital positions required? To solve this, let n be the 
base and m the number of places (decimal or otherwise) required. The 
minimum value of D = nm subject to the condition n™ = K has to be 
found. An easy application of the calculus gives the value 


it é. 


€ 


The use of this transcendental number as a base is impracticable and 2, 3, 
or 4 are logically the next best choice. A numerical comparison is interest- 
ing; to record numbers up to 10° in decimal notation 60 digital positions 
are required, in binary or quaternary scales 40 digital positions, and in 
ternary scale 38 positions. The decimal notation is thus inefficient by a 
factor of about 1-5. 

It is, of course, possible to adopt a ‘binary decimal’ scale, in which 
decimal digits are represented by their binary equivalents. Thus 8 would 
appear as (1000), 68 as (0110), (1000) etc. Using this representation, 48 
positions would be required to represent numbers up to 106 and the scale 
is therefore inefficient by a factor of 1-2 compared with the binary notation. 

The choice between binary and ternary scale is dictated by the fact that 
the binary representation lends itself to easy application of ‘relays’ and 
flip-flops’ (2)—the electronic equivalent of the relay. 

Une objection frequently raised against binary scale is the difficulty 
experienced by humans in the reading of large numbers of ones and zeros. 
As will be shown later, however, a properly designed machine can take its 
input in decimal form, convert to binary for its own operations, and then 
reconvert its output to decimal form. 

Another point must be mentioned in this connexion. Binary addition 
consists of three possible operations (1+1 = 10, 1+0 = 1, 0+0 = 0), 
whereas in decimal notation there are fifty-five different operations. 
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Multiplication is correspondingly more complicated. The engineering prob- 
lems to be solved in constructing a decimal machine are therefore much 
more formidable than those encountered when using binary scale. This 
objection applies equally to the ‘binary-decimal’ scale. 


Serial and parallel operation 

The next point which arises in the design of a computer is the form of 
operation—serial or parallel. In the first type the digits of any number 
become available one at a time starting from the least significant, whilst, 
in the latter, all digits are simultaneously available. The decision as to which 
form of operation is to be used in any particular machine is one which can 
be made only when the form of the ‘memory’ of the machine is known, and 
leads naturally to a discussion of the latter. 


General properties of the memory 

It was mentioned earlier that human agencies are too slow to control 
machines which are much faster than those already in use. Some means 
of ordering the operations of the machine at speeds greatly in excess 
of neural reaction times is therefore required. This suggests that the 
machine should have a ‘memory’, and the characteristics of this will now 
be considered. 

In the past, machines have been built in which the memory consists of 
two distinct parts, one for orders and the other for numbers. This is. 
as will be shown, an unsound arrangement, both from the viewpoint of 
logical arithmetic and from a much more elementary aspect. Two general 
classes of problem exist: those in which a small number of orders suffice 
for the handling of a large amount of numerical data, and those in which 
a large number of orders operate on a few numbers and produce a single 
number as a solution. If only a finite memory space is available, any idea 
of dividing it into two independent halves for orders and numbers respec- 
tively and exclusively is untenable on these grounds alone, and it will 
therefore be assumed that both orders and numbers are to be placed in 
the memory and a selection mechanism provided to pick out the correct 
material. This leads to the idea of a code. 


Basic components of a computer 

If orders are to be stored in the same memory as numbers, they must 
be put into numerical form. This is relatively simple, but requires a careful 
consideration of the precise set of orders on which the machine is to operate. 
Before considering this, however, the operating components of the machine 
must be defined. Fig. 2 gives a general view of the whole. 
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The main sequence of operations is as follows: 

(1) Control extracts order from memory. 

(2) Control decodes order and directs memory and arithmetic unit to 
execute it. 

(3) Arithmetic unit emits signal that its operation is complete. 

(4) Control locates next order in sequence and extracts it from the 


memory. 





Fic. 2. 


Subsidiary sequences allow the control to jump to other orders not in 
sequence; these shifts may be either conditional or unconditional (vide 
i nfra ). 

The exact constitution of an order is of some importance; it is found that 
a satisfactory make up is: 

‘Perform operation A using data in memory position B.’ 

It is not necessary to specify the place of disposal of the result or the memory 
location of the next order. It will generally be convenient for the latter 
to be in sequence for the greater part of a calculation so that an automatic 
advance of the control (via a counter) performs this sequencing operation. 
For out-of-sequence orders special instructions are needed: 

‘Control execute order contained in memory position K and then 

proceed in sequence from this order’, 
which is the unconditional transfer referred to above; or 

‘If x <0 execute order contained in memory position A and then 

proceed in sequence from this order’, 


which is the conditional transfer. 


The serial memory 

More detailed description of the memory and arithmetic units depends 
almost entirely upon the details of the former. Serial memories suggested 
up to the present depend almost universally on the delay-line principle, 
Fig. 3. 

A tube of liquid (usually mercury) has at its ends two quartz crystals 
, and Y,. On applying an electric impulse to Q,, its piezo-electric property 
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causes it to emit a mechanical pulse which is propagated down the mercury 
column as a compression wave. This compression wave affects a second 











Q, Q, 
Ith Hg rit] 
Fria. 3. 


quartz crystal @, which this time emits an electrical impulse; this impulse 
is led back to Q, (after suitable amplification and reshaping) and the whole 
cycle repeats itself. 





Bp—A—ArAs| {Gro | ass 











Carry 
Delay 
Fig. 4. 


The tubes in common use have a delay of about 1 millisecond so that if 
pulses are applied to Q; at megacycle rate, up to 1,000 can be stored in 
the tube. A binary number now consists of an electrical impulse pattern 
of the type 


The disadvantage of the scheme is of course that, on the average, } milli- 
second will elapse before any given number in the tube becomes available. 
To perform addition the two numbers are taken from different delay 
tubes and combined to form the sum, which is recorded in one of the original 
tubes. The addition of two binary digits may produce a carry; this has to 
be delayed by a time-interval A and fed into the next pair of digits to 
become available, A being the interval between pulses (see Fig. 4). 

By a slightly more complicated arrangement, involving the storing of B 
in a short delay line, multiplication and division can be performed. 

Professor F. C. Williams of Manchester has also developed a storage 
device. He uses a standard cathode-ray tube and stores data as charged 
areas on the screen, zero and one being represented by dots and dashes. 
Numbers are read off by scanning the rows of information, and digits are 
therefore obtained serially. Professor Williams has succeeded in storing 
2048 digits in a 12 in. tube and in retaining information for what 1s 
effectively an indefinite period with no signs of deterioration. 
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The parallel operation memory 

Of the memories which have been or are being developed for parallel 
operation machines one of the most promising is the R.C.A. ‘Selectron’. 
In essence this retains data by storage of charge on a dielectric medium. 
The mode of selection of the particular region for a given digit is interesting. 
A filament (Fig. 5) emits electrons in the direction of the dielectric medium 
which is shielded by two parallel grids of mutually perpendicular wires. 


(a —- + = 
a Transparent 
“to electrons 
r 
4 a 
% + 
~ 
= . 
3 fe 
Fia. 5. 


Normally the grids are negative with respect to the filament; to select a 
particular region of the dielectric two adjacent wires in each grid are made 
positive; the region so defined becomes transparent to electrons, which can 
thus pass through and affect the dielectric storage medium. 


The arithmetic unit for parallel operation 
In the parallel operation machine the digits become available at the 
same time, and different arrangements for addition, multiplication, and 
division have to be adopted. It is generally agreed that the arithmetic unit 
must consist of: 
(1) An accumulator; that is, a device which willadd an incoming number 
to that already in it and store the result. 
(2) A register to store the multiplier in multiplication and the quotient 


in division. 




















In multiplication and division the multiplicand and the divisor are stored 
In a register associated with the output from the memory. Fig. 6 gives 
a diagrammatic representation of the kind of arithmetic unit envisaged. 
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The register and accumulator have facilities which enable numbers con- 
tained therein to be shifted to the right or left. Interconnexions, indicated 
by full lines, are provided between the memory register (/.R.), the shifting 
register (R), and the accumulator (A) to enable numbers to be transferred 
from one unit to another. Numbers may also be transferred from the register 
and accumulator directly to the memory. 

The multiplication sequence (positive numbers only) is as follows: 

If right-hand digit of the multiplier in the register R is unity add M.R. 
into A and shift contents of A one place to right, the right-hand digit of 4 
going to the left-hand digit position of R whose contents also shift one place 
to the right (see dotted lines in Fig. 6). The right-hand digit of R is lost 
and the next digit takes its place as the multiplier for the next cycle. 

If the right-hand digit of R is zero the contents of A and R are simply 
shifted as before but without addition from M.R. 


The accumulator 

It is interesting to examine the mode of operation of the accumulator. 
Just as in the serial type machine, the adding unit has two inputs together 
with a carry input from previous stages. 





This carry has been a source of speed limitation in many projects, since 
the speed of operation of the circuit has to be adjusted so that a carry (in 
the worst possible case) originating at A, B, can be propagated through the 
whole chain of adding units. This limitation is not, however, essential and 
the following circuit remedies the defect: 


Bn 





Here, if an incoming carry pulse to any stage (say A, B,) is going to 
produce a carry (i.e. if the result of the (A; B,) addition is 1), the gate @; 
is closed and the incoming pulse by-passes (A, B,) and proceeds unhindered 
along the chain until it meets a stage which is not going to produce a further 
carry. 
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con- The code 

ated Enough has now been said to make clear the general mode of operation 
ting of the machine; it is proposed, however, before describing individual 
Tred machines, to consider the question of the code in more detail. The following 
ister remarks will apply more to the parallel operation type of machine as it 


seems likely to replace the serial variety as being faster by at least an 
order of magnitude, and it is significant that Professor von Neumann, who 


f.R. pioneered the serial machine, has forsaken it in favour of a parallel operation 
of A type. 

lace If the idea of programming all arithmetic out of the simplest operations 
lost is accepted, the list of orders will contain at least: 


(1) Left shift. 
nply (2) Right shift. 
3) Control transfer. 


(4) Memory to accumulator. 
(5) Accumulator to memory. 


hor (6) Control transfer to order at memory location (2). 


ther Out of this simple set all arithmetic operations can be generated, although 
sreater detail will not be given here save to mention, as an example, that 
the nullity of a number can be sensed and a conditional control transfer 
generated from this simple set. This is simply done by successively 
separating (by suitable shifts) the digits, of the numbers whose nullity is 
required, and substituting into an order which sends the control to memory 
location (x); if the digit being substituted is different from zero, (x) is 
increased by unity and the order found at location (x+-1) forms the start 
= ofa new sequence. In practice a more sophisticated code is adopted and 
7 (in 
the 
and 


1s an example that of our own computer A.R.C. is given. The code for 
iserial type machine may be different in form and for details reference may 
be made to the report on A.C.E. (3). 

Although it has been shown how conditional transfer could be generated 
trom logically simple operations, it is a great convenience to have this 
operation available in the code. It is interesting to consider its best formula- 
tion; two possibilities spring to mind: 

(1) Transfer if X < 0. 

2) Transfer if X 0. 





x to Une can be generated from the other, thus using form (1) consideration of 
> Gs ~ X| gives (2), whereas given (2) a consideration of |X|+X gives (1). 
red From the engineering point of view (1) is more suited to a machine which 
her represents negative numbers by their complements. The question of the 
representation of negatives is an important one. In some machines where 
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magnitudes are restricted to |X| <1 it) turns out that representation ] Aiker 
by complements is a valuable addition since it allows the use of the The 
number —1. Aiken 
4 " ~ ‘ 
Cope or A.R.C. puter 
; - . a then 1 
No. | Code Symbol Description h 
—————_—_————_- ine Ra) ME ——— —————. opera 
r | 0,0000 T; to M Fill high-speed memory from input tape. i 
a | 0,0001 T to T; Start tape moving to position 7; and proceed Harv 
with next order. Aik 
3 0,0010 nl; to i+7 to M, to1+7 | Read material on nth tape between i and i The v 
into M. ai 
4 0,001 M, to 1+] to T, Punch material from memory location 1 to ing V 
1+] on to nth tape. al 
: ‘ ee opera 
5 0,0100 Cto M Shift control to order located at M(x). Mt 
6 0,101 Ce to M If A > o shift control as in 5. rh 
7 0,01IO L(N) Shift contents of 4 and R, N places to left. and ¢ 
So that, e.g., 4(0)—A(20); R(o)—R(20) to 
A(o), A(2)—A(20), 0; R(1)—R(20), A(1). mack 
8 | 0,011 R(N) Shift contents of A, N places to right so that, multi 
e.g., A(o), A(1)—A(20) to A(o), A(o), A(1)— 
A(19). ubov 
9 0,1000 M to cd Ai 
10 0,1001 |.M| to cA ace 
II 0,I010 —M tocA =p 
12 0,101 —|M|tocA unit 
13 0,1100 M to A will. 
14 0,IIOI ‘ |[M|to A 
15 0,III0 —MtoA to Ui 
16 O,1III —|M|to A ; 7 tane 
17 | T,0000 M.R to cA Clear A. Multiply / by R, place ist 20 digits 
| and sign in A after adding unity to 2zst the | 
digit. Leave last 20 digits in R. 5 tinu 
18 1,0001 M.R to cA(N.R.) Perform multiplication without rounding off. E 
19 | I,0010 A-~-MtocR Divide A by M, place quotient in R with last ? 
digit unity. Leave remainder in A. shov 
20} I,00II M tocR 
(e.g. 
21 I,0100 R to cA : 
aa I,O101 RtoM | seen 
23 I,O110 Ato M ie ; mat 
24 I,O1r! A;,toM Substitute digits 1-8 of A in order located 
at M(x). and 
25 I,1000 AtocR T 
26 I,1001 E Signal completion of operation. 
the 
An. 
Computer projects in America pro 
We shall now give a short survey of current calculating machine projects scal 
in America. No description will be given of the E.N.I.A.C., as Hartree has inte 
provided this in his book Calculating Machines (C.U.P.). In any case, of 7} 
although this machine is of historic interest as the first large electronic and 
‘alculator to operate successfully, it has been rapidly outmoded by recent 1 
developments, and cannot be regarded as typical of modern ideas on sior 
automatic computers. of ¢ 
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Aiken’s relay computers 


The A.S.C.C. (Automatic Sequence Controlled Calculator), built by 
\iken at Harvard, also falls into this category. This relay controlled com- 
puter was completed in 1941 and has operated almost continuously since 


then to the present day. A very full description of this machine and its 


eration is given in A Manual of Operation for the A.S.C.C. published by 
Harvard. 

Aiken has now completed a second, and larger, edition of this machine. 
[he work of construction and wiring was completed last June (1947), test- 
ing was in progress, and the arithmetic unit had just been put into 
yperating order. 

The machine consists of the control and four arithmetic units. Input 
ind output is via punched tape and results can be printed directly by the 
machine. It performs the elementary operations of addition, subtraction, 
multiplication, and division, and has the discrimination order described 
ibove. 

Aiken hopes to attain a multiplication time of the order of 1 sec., an 
improvement by a factor of 5 over the A.S.C.C., and a multiple arithmetic 
init has been incorporated with the idea of speeding up the machine. This 
will, of course, greatly increase the complexity of programming, as in order 
to use the machine to best advantage all units must be working simul- 
taneously. Indeed, Professor Aiken has mentioned as a major problem 
the difficulty of programming sufficient work to keep the machines con- 
tinually busy 

Experience of programming for our own relay machine, A.R.C., has 
shown, however, that, using the code given above, quite complex problems 
e.g. matrix multiplication) can be programmed in about 30 minutes, and it 
seems reasonable to suppose that most calculations could be dealt with in a 
matter of hours. It may be concluded from this that an increase in the speed 
ind simplicity of the arithmetic unit is preferable to an increase in size. 


The chief reason for the use of binary scale for automatic computers is 
the lack of a reasonably compact and economical form of decimal memory 
\n essential difference between Aiken’s relay machines and other current 
projects is that they have only a very small memory and operate in decimal 
scale. Orders are taken in sequence from punched tape, and all results, 
intermediate or final, are recorded on tape. This increases the complexity 
of programming and reduces the speed of the machine, as reading from 
ind recording on tape are comparatively slow operations. 

The use of decimal scale, of course, eliminates the necessity for conver- 
sions at the beginning and end of every calculation and reduces the number 


of digits which must be carried, but it has the disadvantage that, from the 


5092.6 
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engineering point of view, the operations of arithmetic are much more 
complicated in decimal than in binary scale. 

In spite of these defects, however, both relay machines are fine examples 
of engineering, and Mark I has a most impressive record of continuoys 
performance. 

Aiken is also building an electronic computer, but development was not 
very far advanced last summer (1947). The main outlines of the machine had 
been decided upon, but the detailed circuits were not then designed. This 
machine is to have a high-speed memory, consisting of a drum with a 
number of magnetic tracks, each with its recording and pick-up head. 
Here again, decimal scale is to be used, decimal numbers being represented 
by a binary code. 

The Bell relay machine 

The only other large relay computer in America has been constructed 
by Bell Telephones, and is now at the Aberdeen Proving Ground, Mary- 
land. This machine is a general purpose calculator, built to the basic plan 
indicated earlier in this report. It uses modified binary scale, has tape 
input and output, and uses a code similar to that outlined above. It has 
been successfully used to solve various differential equations, and for the 
calculation of a table of the binomial function. The limitations of this 
machine lie in the size of the high-speed memory—twenty numbers, stored 





in relay banks—and the experience of those who have used it shows that 
this restriction in memory size both reduces the effective speed of the 
machine and increases the complexity of programming. 

As is projected for other machines using binary scale, the Bell machine 
performs its own conversions to and from this scale. It has a certain degree 
of internal checking in that, when information is transferred from one part 
of the machine to another, the number of 1’s occurring in the number 
before and after transfer are compared; if these are found to be different, 
the machine is automatically stopped. This is by no means a complete 
check, and indeed, no current project has yet found a really satisfactory 
solution to this problem. It has been proposed to run two identical 
machines, coupled in parallel, and so arranged that if the results produced 
at any stage are not identical, the calculation is stopped. This, however, 
is likely to be costly in practice, and, in any case, is not an absolute check 
since the comparing mechanism is liable to error. 

The. E.D.V.A.C. 

So far only ‘parallel operation’ machines have been described. Next 
comes the only overt projectt in America building a serial type compute?, 

+ Eckert and Mauchly’s ‘‘Univac”’ is a commercial project and no details are available 
for publication. 
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the E.D.V.A.C. (Electronic Discrete Variable Automatic Calculator) at the 
Moore School of Electrical Engineering, Philadelphia. This is to be an 
all-purpose, electronic computer using as memory, mercury delay-lines. 
Its counterparts in England are Wilkes’s E.D.S.A.C. at Cambridge and 
the National Physical Laboratory’s A.C.E. 

In July (1947) the main body of E.D.V.A.C. was still under construction. 
\ good deal of work has been done on delay-lines, and in particular on the 
problems arising from corrosion of the tanks by the mercury, and a satis- 
factory model has been adopted. The designs for the accumulator and 
control are complete, and a pilot model of the arithmetic unit was operating. 
This model, ‘Shadrac’ by name, performs addition, subtraction, multiplica- 
tion, and division, but no controls were built in and numbers had to be 
inserted by hand. ‘Shadrac’ was built for experimental purposes only and 
has a capacity for twenty digit binary numbers, the delay lines being about 
I8in. long. The final machine will work to forty binary digits and on the 
basis of experience gained with ‘Shadrac’, it is hoped to achieve a multiplica- 
tion time of 1 millisecond. 

A typical order will consist of four parts: 

(1) and (2) Specification of the locations of the two numbers to be 

operated upon. 

(3) Operation (+, etc.). 

(4) Location to which result is to be sent. 

Some coded examples are given in a report on E.D.V.A.C. (not published) 
ind these indicate that coding for the machine is a comparatively quick 
and easy process. Orders and data will be punched on teletype tape, 
transferred to magnetized wire, and thence to the machine. This rather 
complicated process is necessary since direct tape input would be far too 
slow for the machine. The same method is to be adopted for the project 
next to be described, that of von Neumann at the Institute for Advanced 
Study, Princeton. 


The Princeton electronic computer 

This is probably the most ambitious project at present in existence, as 
it is hoped to attain a multiplication time of less than 100 sec. The 
machine is to be a parallel type electronic computer using binary scale. 
[t differs from all those discussed so far in the size of its memory; this is 
hoped to be about 4,000, forty digit binary numbers to be stored in the 
‘selectrons’ briefly described above. 

The machine will follow the general plan described earlier in this report, 
and will consist of the memory, accumulator, shifting register, control, and 
the input and output units. 
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The input is to be via teletype tape and magnetized wire, and the 
machine will be provided with a library of these wires containing tabulated 
functions and codes for routine calculations. 

The code is again similar to that described above, and each order will 
consist of two parts: 

(1) The operation to be performed. 

(2) The memory location of the number to be operated upon. 
Conversions to and from binary scale will be done by the machine, decimal] 
numbers being inserted in binary coded form. The large memory capacity 
and simple code, will make the coding of problems for this machine 
comparatively rapid and easy. 

At the end of last summer a shifting register and accumulator were in 
operation, and the magnetic input and output was complete. The accumu- 
lator had an addition time of 2 sec., but it was hoped to improve on this 
the control circuits necessary to make multiplication and division possibk 
had not, at that time, been completed. 

A great deal of work has been done on the best and most economical 
medium for the magnetic input and wire has been chosen as providing 
minimum bulk with good recording capacity. To avoid undue bulk, small 
gauge wire will be used, and difficulty was experienced at first because the 
tensions set up when starting and stopping caused breakages. This problem 
has been solved by mounting two drums on a common shaft and winding 
the wire from one to the other via the reading and recording heads. The 
effective radii of these drums will, of course, vary slightly with the amount 
of wire on them, but the consequent tensions in the wire are eliminated 
by the use of a differential gear and a small servo-motor. 

Another problem which has to be considered in this connexion is that 
of reading from the wire while accelerating. It is possible to lose digits in 
this process, and this is to be avoided by having a set of marker pulses 
between each ‘word’; A counting device records the number of digits read 
out between each marker pulse, and if any are missed, the word is rejected. 
It is hoped, in the final machine, to record 100 digits to the inch and to run 
the wire at 50 ft. per sec., giving an input speed of 60,000 digits per sec. 

Progress on the memory has, so far, been slow. Originally it was to 
consist of forty ‘selectrons’, each containing 4,000 digits. At the end ot 
the summer (1947), however, no working model of the selectron had been 
produced, and it seems that, when this is available, it may contain, at first. 
only 256 digits. While this memory is adequate for a relay machine, it will 
probably be quite insufficient for the much faster electronic computer, and it 
may be necessary to develop an alternative form of high-speed memory) 
consisting of, for example, magnetized drums. 
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Forrester, of the Massachusetts Institute of Technology, is engaged on 
an almost identical project to that of von Neumann. Development is 
slightly behind that at Princeton, but the only essential difference is in 
the proposed form of the high-speed memory. Forrester hopes to use some 
form of electrostatic storage, presumably of the general type suggested by 
Williams. 

Althoug] 
computer projects, it is interesting to note that all are sponsored, in part 


the United States is comparatively rich in computers and 


it any rate, by the Government, and that all, with the exception of the 
Princeton machine and Forrester’s computer, are eventually to be moved 
to the Aberdeen proving ground or to the U.S. Navy Proving Ground at 


Dahlgren. Virginia. 
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SUBHARMONICS IN A NON-LINEAR SYSTEM WITH 
UNSYMMETRICAL RESTORING FORCE 
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(Department of Mathematics, The University, Manchester) 


[Received 11 May 1948] 


SUMMARY 

The behaviour of an oscillatory system which is positively damped but in whicl 
the restoring force is unsymmetrical about the equilibrium position is investigated, 
It is shown that if a periodic external force whose period is approximately half th 
natural period is applied, the system may oscillate with twice the period of th 
external force (such oscillations are called subharmonics of order 2). The gsub- 
harmonics only occur when the amplitude of the external force reaches a certain 
critical value ; if the forcing period is slightly greater than half the natural period, 
there is a second critical amplitude beyond which the ordinary forced oscillation 
(with the same period as the external force) becomes unstable. In the range between 
the two critical amplitudes, therefore, the forced oscillation or the subharmonics 
may occur (which of them occurs will depend on the initial displacement and 
velocity of the system), whilst beyond the second critical amplitude only the 
subharmonics can occur. “This gives rise to a ‘hysteresis’ effect when the amplitude 
of the external force is increased and then decreased again. 


1. Introduction 

It is well known that non-linear differential equations of the form 
¥+-f(x)é+-g(x) = f cos pt 

(dots denoting derivatives with respect to ¢) may have periodic solutions 

whose least period is an integral multiple 2nz/p of the period 27/p of the 

forcing term (1). (Such solutions are called subharmonics of order n.) It 

has been usual to assume that the restoring force is symmetrical about 

« = 0, ie. that g(x) is an odd function of x; the simplest equation with 

ain unsymmetrical restoring term is 


+ hat +- w2x(1-+ ax) fos pt. (1) 


This equation (with a positive damping term, i.e. k > 0) arises in the 
theory of loud speakers (2), and has been investigated both experimentally 
and theoretically by Pedersen (3). As, however, Pedersen’s treatment is 
not rigorous and his results appear to be incorrect, we shall discuss the 
existence of subharmonics for equation (1) rigorously by Poincaré’s method 
of small parameters (4). It will be shown that, when & and « are both 
small and p is approximately equal to 2w, equation (1) has under certain 
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conditions subharmonic solutions of order 2. The stability of these solu- 
tions will also be discussed. 
[ should like to thank Miss M. L. Cartwright for suggesting this investi- 


gation to me, and for her advice during its progress. 


2. We may arrange, by changing the time scale, that the forcing term in 
1) is feos 2¢. In practice, w would be held fixed and p varied, but since 
nly small variations of p about p = 2w will be considered they may be 
taken into account by varying w. We assume therefore that p = 2, and 
that w2—1, k, and « are all small and of the same order of magnitude. 
We may then write k = Ax, w* = 1+ pa, and f = 3v, where A, a, and v 
may be assumed to be positive, and A, », and v will be of order unity. 
Equation (1) then becomes 

X+Aaae-+ (1+ px)a(1- aa) 3v cos 2¢. (2) 
It will appear that (2) has subharmonics whose amplitude is of order « 
and it is therefore convenient to write « = e? and y = ez, so that 

ij +-Ae*y+- (1+-pe”)y(1+-ey) 3ve cos 2, (3) 


and the subharmonies will now have amplitudes of order unity. 


3. We shall now prove that, when A, p, and v are fixed, A < v, and « is 
sufficiently small, equation (3) may have subharmonic solutions of order 2 
ie. periodic solutions of period 27). These subharmonics will always occur 
n pairs, since if y = d(t) is a solution of (3), so is y = ¢(t+-7). 

When e = 0, the solution of (3) with initial values y(0) = a, 4(0) = b 
is y = o,(t) = acost+bsint, and has period 27. We shall attempt to find 
solutions of (3), of period 27, which reduce to such a ‘generating solution’ 
when € > 0. 

The solution of (3) with y(0) a, y(0) b is analytic in ¢ for sufficiently 
small e«, and is of the form 

y do(t)+ ed, (t)+-€7,(t)+..., (4) 
where the functions ¢,(¢) are analytic in ¢, a, b, and the series converges 
uniformly when ¢, a, 6 lie in fixed finite intervals and « is sufficiently 
small (5). Substituting from (4) into (3), we find that 

¢,+¢4, 3v cos 2t— gi, 
$,(0) = ,(0) = 0; 


bo+$z Abo—Hbo—2bod1 = yf, say, | 


(6 
$(0) = (0) = 0. J 


(5) 


From (5) we obtain 
$,(t) }(a?+-b?)+-(v+-4a?-+ 2b?)cos t— Zab sin t— 


—(v—}a?+ 4b?)cos 2t+ hab sin 2¢, 
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and hence we may calculate s(t), the right-hand member of (6). We shal 
only require the values of ¢,(27) and ¢,(27), and from (6) these are given by 
on 


$,(277) [ b(t)sin t dt, 


on 
ys(t)cos t dt. 

0 

Inserting the value of y(t) obtained from (5) and (6) (only terms involving 
cost or sint being required) we finally obtain 


Vi 


Lg (2m) Aa+-(p+v)b—sb(a?+-b?) = P(a,b), say, | 


l 55(2rr) Ab—(p—v)a-+ia(a*+-b?) = Q(a,b), say. 


Hence, from (4), we find that 
y(27)—y(0) = ze*| P(a, b)+-eP,(a, b)+...], 
y(27)—y(0) = we Q(a, b)+-€Q, (a, b)+...], J 
where P(a,b), Q(a,b) are given by (7), and P,(a,b), Q,(a,b),... are poly- 
nomials in a,b (whose values will not be needed further). Our task is now 
to choose a and b (as functions of €) so that y(27) = y(0) and 4(27) = ¥(0). 
To this end, we first choose a, and b, so that 


P (ao, 59) Q(4, 59) = 0 


(8 


and then write a = a,+&, b = b)+ 7. By Poincaré’s expansion theorem 
(7), y(27)—y(0) and 4(27)—y(0) will then be analytic functions of &, », « 
which can be expanded in power series convergent when |€|, |7|, and « 
are sufficiently small. These expansions can be obtained by expanding 
P(a,b), P,(a,b),..., Q(a,b), Q,(a, 5),... in powers of € and », and inserting 
these expansions in (8); rearranging in powers of € and », we obtain 
expansions 


ane an A(e)+ Ble)E+Cle)n+... = X(E, n, €), say, 


9 
we 





(9) 
a 
y2n) —W\0) = D(e)+ E(e)E+ F(e)n+... = V(E, n, €), say, 
Tre“ 
where A(e),..., F(e) are analytic for sufficiently small «. Using (7) and (3), 


we find that 
A(0) = D(0) = 0, 
B(0) —A—$ayby, C(0) ptv—saz—Sb5, 


E(0) = —p+v+$ag+9bp, — -F(0) = —A+ jag bo. 
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a(x, Y) 
Olé, ”) 


_ | B(O) C(0)| 
| E(0) F(0)| 


Hence : A, say, 





é n= 0 
snd so, provided A + 0, the equations 
X(é, nN; €) Y (€, n, €) — () 


have, for sufficiently small ¢, a solution (analytic in e) 


s E(e), n nile), 
with (0) = »(0) = 0. The solution of (3) with 
y(0) a,+€(e), y(O) bo + ne) 


then has period 27 and is of the form y = y(t,e) (analytic in ¢), with 


y(t. 0) a, cos t- b,sint. 


4, The equations P (ao, bo) Q(d9, 69) = 0 


which determine the possible subharmonics will now be discussed. The 
solution dy = b 0 corresponds to a periodic solution of period 7, which 


0 


will be discussed in § 5. Any other solution must satisfy 
tlt = tA 2 a0 
+05 = Murer}. J 


These equations have solutions (other than a, = b, = 0) if and only if 
\<v and p a|(v?—A?*). If -4/(v?—A?) <p a/(v?—A?) there are 
two solutions, with the lower sign in (10), and if p > J(v?—A?) there are 
four solutions, with both signs in (10). Each sign in (10) yields a pair of 
solutions, (@,b)) and (—ay, —b,), corresponding to a pair of subharmonics 
with phase difference 7z. 

With the values of a, and by given by (10), it is found that 

A 4,/(A2—v?){,/(v?—A?) Fp}. 
Since the upper sign occurs only if » > ./(v?—A?*), and the lower sign only 
If pe J(A2—v?), A 0 with either sign provided A < v, and this finally 
establishes the existence of subharmonics for equation (3). 

Equation (9) also determines the stability properties of the sub- 
harmonics. If a subharmonic has initial values y(0) = a, y(0) = 6, then 
the solution with initial values a+€ and b+7 has y(27) = a+’, 
(27) = b+-»', where from (9) 

é’ = [1+ we? B(0) |€+-7e2C(0)n+..., 
7’ me" H(O)E4 [1 L me*F(0)|n-4 ver 


the terms omitted involving second and higher powers of £ and » or third 
“« 
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and higher powers of «. The stability of the subharmonic therefore 
depends (7) on the characteristic roots p of the matrix 
1+ ze? B(O) meO(0) \. 
me? E(0) 1+ ze? F(0))’ 
these roots are p 1+-e*{—A-+,/(A?—A)} (using the values of B(0),..., F(0) 
found in § 3). 
There are three cases to discuss: 
(a) The periodic solution (ay = by = 0) has A = A?+-p?—v*, 
p 1+-e*{—A+,/(v?—p?)}. 
If A >v or if \<v and |p| > J(v?—A*) both roots have modulus less 
than 1 and the solution is completely stable; if A < v and |u| < ./(v?—)’) 
both roots are real, one being less than 1 and one greater than 1, so that 
the solution is ‘directly unstable’ (of ‘col’ type). 
(b) A <v, p > —,(v?—A?), lower sign in (10). Both roots have modulus 
less than 1 and the solution is completely stable. 
(c)A<vp> J(v? —A*), upper sign in (10). Both roots are real, one 
less than 1 and one greater than 1; the solution is directly unstable. 
5. It was stated in § 4 that the solution a, = b, = 0 of 
P(d 9,69) = Q(ao; b9) = 0 
corresponds to a periodic solution of period 7 of equation (3); the ampli- 
tude of this solution will be of order a?. To prove these statements, a 
method similar to that of § 3 may be used, working with equation (2) and 
expanding in powers of «. It is easily shown that the periodic solution 
exists for small «, and is of the form 
x = a(t,a) (analytic in «) 


with x(t, 0) vcos 2t. We shall omit a detailed proof. 


6. Weshall now summarize the results of §§ 3 and 4, applied to equation (2): 
¢+kei+w'r(1+- ax) = 3vcos 2¢. 
The results are valid when « < ap, where «, depends on k/a, (w?—1)/a, 
and v; in practice this will mean that «, k, and w?—1 must be small and of 
the same order of magnitude and that v must not be too large. 
(a) There is always a periodic solution of the form 
x —v cos 2t+O(a); 

this solution is completely stable unless k < va and |w?—1| < (v?x®—h*). 
when it is directly unstable. 


(b) There are subharmonies if k < va and w? > 1—,/(v?a?—k’). There 
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are two or four subharmonics, occurring in pairs with phase-difference 7, 


according as |w*— 1} (v2x?— k?)? or w? > 1+,/(v? x? —k?) respectively. 
(c) When k < va and |w?—1| < ,/(v?a?—k?) the two subharmonics are 


completely stable, and are of the form 


x p cos(t 6)+-A cost+ Bsint 


> 


}ap?—v cos 2t— dap? cos(2t+ 26)+-O(a?), (11) 


> v.99 9 
: 6 9 .* " va-t-./(v?«?— k? 
p* {w?— 1-+-,/(v?a?— k?)} and tan @ esi rally" 





giving two possible values of @ differing by 7). The coefficients A and B 
are O(1) but cannot be determined from our calculations; they are, how- 
ever, usually of smaller order than p, which is of order a~? unless w? is 
very near to 1—,/(v?a®—k?). 

(4) When i vx and w? 1+-,/(v2a?—k?), the stable pair of sub- 
harmonics is still present, but there is also a directly unstable pair, again 
given by (11) but with 


. 6 . 7" ‘ Va (v2oq?—k? 
p* —{w?—1—,/(v?a?—k*)} and tané ¥ E as 
7. It is of some interest to examine how the solutions vary as the ampli- 
tude or the frequency of the forcing term is varied. 

(1) Varying amplitude (k, «, w fixed; v varied): 

(2) As v increases from 0 to k/a, there is a stable periodic solution of 
amplitude approximately v, and there are no subharmonics. 

(b) If w > 1, two pairs of subharmonics appear when v reaches k/a—a 
stable pair with amplitude approximately 
Vox} w?—1- | (v2.0? k?) |}, 


ul 


dan unstable pair with amplitude approximately 


é | 


v8 ¥ Uw? l a (v7? k?) |}; 
the periodic solution remains stable. 

When v reaches a~1,/{k?-+-(w?—1)*}, the unstable pair of subharmonics 
disappears, the stable pair remains, and the periodic solution becomes 
unstable. 

: . ° ° | . . ‘ 

(c) Ifw < 1, no subharmonics appear until v reaches a-1,/{k?+ (w?— 1)*}; 
a stable pair of subharmonics, with amplitude approximately 


VE x Vw ] t-4/(v?a?—k?) |, 
then appears and the periodic solution becomes unstable. 
[t will also be seen, from (11), that the subharmonics are approximately 
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sinusoidal when they first appear, but that when v is increased they become 
unsymmetrical about x = 0 and also contain a marked second harmonic. 


(2) Varying frequency (k, «, v fixed; w varied): 


(a) If k >va, there is a stable periodic solution and there are no 
subharmonics. 

(b) If k < va, there are three different ranges of w to consider, If 
w2 1—,/(v?a?—k?), there is a stable periodic solution and there are no 














subharmonics. If |w?—1 V(v?a?—k?) the periodic solution is unstable 
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and there is a stable pair of subharmonics with approximate amplitude 
V8 al w?— 1+4-,/(v2a02—k?)]! (but this approximation is invalid when ? is 
very near to 1 —,|(v?a2—k?)). If w? > 1+,/(v?a?— k?) there is a stable 
periodic solution, a stable pair of subharmonics as before, and an unstable 
pair, with approximate amplitude v$a—![w?—1—,/(v?x2 


k?) |! when w? is 
not too near 1-+-,/(v?a?—k?). 

Curves of amplitude p plotted against v (with w constant) or w (with 
v constant) are sketched in Figs. 1 and 2, unstable solutions being indicated 
by dotted lines.+ (It should be noted that these curves are only accurate 


+ In Figs. 1 and 2 s.h. stands for subharmonic and p.s. for periodic solution. 
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when p> 1.) Fig. 2 shows that there is no definite resonance near w = 1, 
but that the amplitude of the subharmonics appears to increase indefinitely 
13 w is increased (i.e. as the frequency of the forcing term is decreased 
below twice the natural frequency of the system). It should, however, be 
borne in mind that our results cease to hold if w?—1 becomes too large 
ompared with «, and that then the first approximation to a subharmonic 
would be a periodic solution of #+w?a(1- xX) 0 rather than of 
Y 0. 

Inspection of Fig. 1 also reveals the possibility (when w > 1) of a ‘hys- 
teresis’ effect if vy is increased and then decreased again. As v increases, the 
stable periodic solution will be followed until v reaches a 1 /{k?+-(w?—1)?}; 
the system will then ‘jump’ to one of the stable subharmonics. If now v 
s decreased, this subharmonic will remain until v reaches k/a, when a 
imp back to the periodic solution will take place. The process is indicated 
vy arrows in Fig. 1. 

A possible objection to using equation (1) for a system with unsym- 
metrical restoring force is that w*a(1+ax) may cease to represent the 


restoring force when |aa| becomes comparable with 1, since the actual 


restoring force probably does not change sign for ax 1. However, 
for the solutions which we have found |aa|<1, and so if necessary the 
restoring term could be modified when ax O(1) without affecting the 


existence of subharmonics. 


8. The results of § 3 also show how any solution of equation (3) approaches 
1 subharmonic or periodic solution. For this purpose, we return to 
equation (8), which shows that the solution with y(0) = a, 4(0) = b has 
(27) a+t-Aa,. 7(27) b-+-Ab, where 

Aa/zre* P(a,b)+- Ole), | (12) 

Ab/me2 = Q(a,b)+ Ole). J - 
The behaviour of solutions of (3) as ¢ increases by multiples of 27 is then 
ndicated by the behaviour of points of the (a, b)-plane under iterations of 
the transformation 7’: (a,b) > (a+Aa, b+Ab). Approach to a sub- 
harmonic or periodic solution corresponds to approach of (a,b) to a point 
M,b9) with P(a9,b5) = Q(ao, by) = 9, i.e. to a fixed point of the trans- 
lormation 7’. By drawing the curves P(a,b) 0 and Q(a,b) 0, the 
4,b)-plane can be divided into regions in which Aa 2 0 and Ab 2 0, and 
lence the shape of the ‘trajectories’ described by points (a,b) under 
iterations of 7’ can be estimated. Two typical configurations are sketched 
in Figs. 3 and 4; and signs of Aa and Ab are indicated by symbols like 
| (to denote a region in which Aa 0 and Ab > 0). 
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In Fig. 3 there are two stable fixed points S,, 8, of ‘focus’ type, and one 
unstable fixed point C of ‘col’ type. It is known (8) that there are two 









ee wee 
noe 


- 
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curves invariant under 7’ passing through a fixed point C of ‘col’ type: 
points on one of these curves approach C and points on the other recede 
Ms ¥ . . ~ 7 . . 4 

from C under iterations of 7’. These curves are indicated by broken lines; 
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the unstable curve spirals into the two stable fixed points, whilst the stable 
urve spirals inwards for large values of a?+-b? before settling down into 
asteady approach to C. 

In Fig. 4 there are three stable fixed points Sp, S,, S, (Sy corresponding to 
the periodic solution) of ‘focus’ type, and two unstable fixed points C,, C, of 
col’ type. In this case, the unstable invariant curve through C,(C,) spirals 
into S, and S,(S,). 

In both cases, points at some distance from the origin spiral inwards 
until they are attracted towards one of the stable fixed points. Which 
fixed point is finally reached depends on the initial position of (a,b); the 
invariant curves through the cols divide any ray through the origin into 
1 sequence of intervals from which points reach S,, S,, S,,... in Fig. 3 and 
§,. S,. Sp, So, So, Sy,... in Fig. 4. Points near one of the invariant curves 
will spend some time near one of the cols (first approaching a col and then 
receding again); the corresponding solution will for a limited time give the 
mpression of settling down to a steady periodic state (periodic solution 
in Fig. 3, subharmonic in Fig. 4) before settling down to a different final 
state (subharmonic in Fig. 3, periodic solution or subharmonic in Fig. 4). 
The parameters in Figs. 3 and 4 have been chosen so that the stable fixed 
points are all of ‘focus’ type: for certain ranges of the parameters, it is 
possible to have fixed points of ‘node’ type, for which there is a steady 
one-dimensional’ approach to the fixed point. The approach of a solution 
to a subharmonic in the two cases is of a different nature (9); for a focus, 
the approach is like p(¢)cos|¢—6@(t)|, where p(t) and @(¢) oscillate about their 
final values; for a node, p(t) and @(t) tend steadily in one direction towards 
their final values. From the results of § 4, it follows that the stable 
subharmoniecs are of ‘node’ or ‘focus’ type according as w? is less than or 


5k? — 4202 — , a 
greater than 14 _.. and that the periodic solution (when it is 
4, /(v?? — k?) 
stable) is of ‘node’ or ‘focus’ type according as |w?—1| < va or 


y" ] VX 
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A NOTE ON THE NUMERICAL INTEGRATION 
OF DIFFERENTIAL EQUATIONS 


By E. M. WILSON 


(Royal Naval Scientific Service, Admiralty, London, S.W.1) 
[Received 19 February 1948] 


THE Taylor-series method for step-by-step integration of a differential 
equation is well known and used, particularly when high accuracy is 
needed. However, attention does not seem to have been drawn to the 
following version of the method, which, though not of universal application, 
can often be applied to linear equations. For an equation of order m, this 
form of the method consists in expressing the solution and the first (m—1) 
derivatives at a point of tabulation as sums of the values of the same 
functions at the preceding tabular point multiplied by previously computed 
functions of the independent variable. 
Consider, first, the linear first-order equation 

rd P,(x)y+Q,(2). (1) 
By repeated differentiation of (and substitution from) equation (1), we 


have 
d’y 
dat =e 


P(x)y+Q,(2), (2) 


where the successive P’s and Q’s are given by the recurrence relations 


LP. 
P= 5° +hP, 
Ga 
K (3) 
dQ. 
Qt “4 QP, 
da 
We may thus write Taylor’s series as 
tw) = (14+0P,4" P4 + (we, + 5 Py(x)+Q 
y(a+w) WE +5 Lot Jy(e)+( WG +5 i+ -- ya 
(4) 


P and Q are functions of x only (for a fixed interval of tabulation), and 
can be calculated before the integration begins. Their calculation may or 
may not be easy, but once they have been obtained, the solution of the 
differential equation can be built up in the most straightforward manner. 
[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 2 (1949)] 
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Correspondingly, for a second-order equation 
d*y ' da ” 
* = P(x)y+Q(x)—2 + R(x), (5) 
dx? dx 
we should have two finite difference equations: 
y(a+w) = P,y(x)+Q,y'(x)+R,(2) (6) 
y (x«+w) = P,y(x)+Q.y'(x)+R,(2) 


where dashes denote derivatives ; and so on for equations of higher orders. 

The number of elementary operations (multiplications and additions) 
required after the calculation of the functions P, Q, R may be sufficient 
to make the process uneconomic, unless a large part of the work can be 
performed mechanically. Until recently, there was no machine suitable 
for carrying out a continuous series of multiplications and additions, but 
the new automatic computing machines would seem perfectly suited to 
the process. 

The fact that the accuracy to which P and Q, etc., must be calculated 
has to be determined at the outset entails no more difficulty than is 
normally encountered in deciding how many figures to take initially in y, 
nd this is a decision which depends more upon the degree of stability 
f the equation to be integrated than upon the particular integration 
procedure used. 

It may happen that calculation of the preliminary functions P, Q, etc., 
enough for the process to be useful even without automatic 


machines. As an example, we take the equationt 


y” y' + ary, (7) 
: being a parameter. Here we have 
y P.y+Q,y’, (8) 
vhere ; = 8 
Fiea P,+axQ, Dd 441,51 | (9) 
, > .* 8 . 
Qa = O4+-9,4-P. = Dd 5,41,2 J 
and > > 
A 1) R=0 \ 
Q, = 0) 0; 1J 
the relations between the coefficients a,,, b,, in these finite sums appear 
it once as: 
A415 {oF 1 )Qy g44 T xb, « 1 (10) 
b, 1,s (s+] )b,.s +177 Ras Tas 


equation can be reduced to that of the Airy integral, Y” = XY. 


P 


5092.6 
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from which we can readily construct the following table: 














TABLE I[ 
= Pia ar,s | ie sae br, 
eS 0 l 2 . 4 0 l 2 3 
0 | l 0 
l | 0 l 
2 0 x 1 0 
3 x x 1 x 
4 aX x x? | 2a+1 2a 0 
5 | x x? + 2a? 5a+]1 3a x? 
6 4a7+a Q9a?+a 3a? | 9a+1 6a?+4a 322 0 
| 9a?*+a Llda®+a G9a3+4a2 3a3 | 10a?+14a+1 21a?+5a 6a? 3 
Equations (6) for this case are then 
1 
y(a+w) = y(x){A +A, ar+A, o22?+...34 
+y'(«){By+ B, ox+ B, a?x?+...} uy 
lb aes 64? > et A Mees fs ‘ 
y' (e+w) = y(x){ Ag+ A, or+ Ay o72?+...}+4 
+y'(x){ Bot By ox+ Bg o®a?-+...} 
where, w,, being written for w”/n!, 
Ay = Wp ta(wstwyt...)+a7(4we+t 9w,+...)+... 
, TT : | 
A, = (Wetwst...)+a(4w,+9w,+ 15w,+...)+... 
| 
Ag = (wy+2ws+ 36+...) +a(9w,+...)+... 
A, = (we+3u,+...)+... 
By = (wy twetw st...) +a(2w,+50,+ 9wet l4w,+...)+ p. (12) 
+a?(10w,-+-...) To ee 
B, = (wg+2w,+3w;+...)+a(6wg+21w,+...)+... | 
B, = (w,+3w,+ 6w,+...)+... 
By = (w.-+...) T eee 


The expressions for A}. and B'. differ from those for A, and B, only in 
having w,,_, in place of w,, (w_, = 0, wy = 1). 

The multipliers of y(a) and y’(x) are infinite series in 2, with coefficients 
which are infinite series in a, these series again having infinite series in v 
as coefficients; but, in fact, calculation is easy. Rates of convergence 
depend upon the choice of w, which may be taken so that P,, P,, Q,, and 
Q, can be treated as polynomials in x of low order, and either built up 
at interval w from a constant difference, or, perhaps, interpolated from 
a building-up at a wider interval (when there would be fewer extra figures 
needed). These preliminary computations could be carried out on, say, 
a National accounting machine, or Hollerith tabulator, and the straight- 
forward integration of y then performed with the help of any ordinary 
desk machine. 
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With an interval w of 0-1, equations (12) become: 


Ay 

A, 0:00517 09181 
A, 0-0° 43376-4 
Az 0-08 14+ 

B, 0-10517 O9181 
B, 0-00017 52557 
B, 0-07 876-4 

B, 0-9104 

Ay 0:00517 09181 
A; 0-10517 09181 
Ay 0-00017 52557 
As 0-07 876-4 

Bo 1-10517 09181] 
4 0:00534 61737 
B, = 0-0° 44251 

Bs 0-08 15-4 


We use this interval of w 
decimals for the case a 
values of 


out in Table IT. 


1 ()-()6 


0-1, 


|(-00017 0918la+0-08 57a2-+... 
0-08 346la+... 


0-0? 2a-+-... 


+0-0° 87628a-+... 


0-0° S8a coe 


L()-08 3461a2-+... 


0-00001 


0-07 13 


0-00035 49365a+ 0-07 


5301a-+... 


74379a+... 


la- 


1470? 
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(13) 


0-1 to carry out the integration with ten 
1, taking the initial 


from 


3 0 to x 


y and y’ to be respectively unity and zero. The results are set 











TABLE II 

P =D A,ofx"|Q, = EB, ax") P, = X Apatx"|Q, = X Bratx’ V(x) yy’ (x) 
00/1 I 7091 10517 17944 O51 70953 | 1°10520 64119 | 1°00000 00000 | 000000 00000 
‘I 8801 35470 156 88394 25 98746 I 70919 51 70953 
2 I2 05127 52996 262 06186 31 33382 14 02796 214 03909 
3 17 22244 70522 367 24329 36 68027 | 48 59255 498 67888 
4 22 ¢ 7 88048 472 42822 42 02680 | 118 27293 918 64482 
5 27 56503 8 05574 577 61665 47 37343 237 31507 | 1488 47558 
6 32 7364 3101 82 80859 52 72014 421 50421 2224 45808 
7 37 90798 40628 788 00404 58 06694 688 35230 | 3144 88557 
8 43 07 58155 893 20299 63 41383 | 1057 31306 | 4270 35298 
09 48 251 75682 998 40545 68 76081 1550 02869 | 5624 09498 
I'0| 1°00053 4 10518 93209 1103 61141 | 1:10574 10788 | 1-02190 61274 | 0°07232 37315 

y(a+-O-] y(x)P, (x) + y’(x)Q,(2) ; y (x+0-1) = y(x)P,(x)+y'(x)Q,(2). 


We may compare the final values with those 


102190 61274 and 0-07232 


9r9 


Viv 


14, 


of y(1) and y’(1) found 
by summing the readily obtained ascending series. These are, respectively, 


Many thanks are due to Mr. D. H. Sadler, for helpful comments on the 
method herein described, and to Mr. B. W. Conolly for assistance with 


the example. 








TABLE OF { Io(x)dx FOR x = 0(0-1)20(1)25 


By M. ROTHMAN (Royal Technical College, Glasgow) 


[Received 15 March 1948; revised 3 September 1948] 


THE solution of the problem of a lubricant or viscous fluid flowing through 
two very close horizontal plates was first given by Michell, whose method 
is reproduced by Stanton (1). 

If we consider the plates to be vertical and fixed instead of horizontal 
we are led, by the same analysis, to the equations: 


P : > Pm (1) 
Pm = (@,,8in mz)/ma, (2) 


where p is the pressure at a point of one of the plates and w,,, satisfies the 


equation 29 , 
C Wm | I CWm ] | l | A 0 € 
fa 5 oh DA ee @, (3) 


where € = mx and k = 12pgq/z, this being a Bessel equation with the added 
constant k/m?. 

These equations are the same as those derived by Stanton for the 
horizontal case except for the constant £/m*, and the solution of equation 
(3) will therefore have the same complementary function as the corre- 
sponding equation of Stanton but a slightly different particular integral, 
the solution in this case being 


. v4 
w - A, 2T(0)+B Edt) (5+ Reed } (4) 


where A,,, B 


The total pressure on the plates is then given by P = || pdaxdz, the 


m> &, and m are constants. 


integrals being taken over the plates. The evaluation of P necessitates 
a knowledge of the values of the integrals 


LA) ay and | Ky(2) gy 
j 2 . x 


0 a 





which may also be written in the forms 
. rd 
| A(t) de = | Iy(x) du—[L,(x) |p (5) 
, i 
(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 2 (1949)] 
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TABLE OF A CERTAIN INTEGRAL 


zx 
. 


ie dx = - | Ko(x)de—[Ky(x) hs, (6) 


. . 


x x 


r K, (x) 
and - ( 


30 that the solution of the problem is dependent on tables of the two 
integrals a 
| Ly(w)dx and | Kj(x)dzx. (7) 
0 0 
In the actual solution of the problem it is 

SA», | Lg(max) d(max) 

. 0 
that is wanted, and this is convergent; but it may happen, as was the case 
for a particular experiment, that the integral exceeded the limit of existing 
tables before its product with the constant A,, was ignorable. 

Miiller (2) has computed both the integrals (7) to x = 16, at which point 
the latter is sufficiently near to its limiting value to make further com- 
putation unnecessary. His computation of the former integral is for the 
range x = 0(0-1)5-2(0-2)16, and as it was found convenient to have a 
slightly extended form of this table for a particular experiment, the table 
has been recomputed to x = 20, the results being given to 8 significant 
figures and at intervals of 0-1 up to this value. The values of this integral 
multiplied by e-* for successive integers is also given up to x = 25. 

Several formulae are suitable for the computation of this integral. 

“2 at 
Firstly : I(x) J, (ix) 5 +. ne -e 


24.12.22 | 
and therefore 
: 2 [x 9 (x 3 > x\> 9 a\7 
| I(x)de = =(5)+=(5) + (2) +5 (3 fine (8) 
: 1\2) ' 3\2) * 5(2!j2\2 7(3!)2\2) | 


Secondly, terms in (8) may not be computed individually, but by means of 


—— 2 (s)"" 2n—1 1 [x i. (9) 
n+l (2Qn+1)(n!)2\2 Qn+1 n2\2) ” 


where 7’ is the nth term in (8). 


? 


Thirdly: 


( I(x) dx bora| L(x) {Ip(x) — Lo(x)} —Ly(x){ (x) — L 1(x)}], (10) 


0 


where L,(x) and L_,(x) are modified Struve functions. 
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Fourthly, the asymptotic expansion of (10): 


r 


. = 1.3 hat $3.5 19.3*.5.7 
ey a | | ees ie 
| Zola) da ( ar 3! 


7 4/ (22x) | 1! 8a ) 3! (8a)8 
0 
2.32 12.32.52 
te! gt weOCOS } 
'* 4s 37. 5? 1 13.3 .3*.6 
14 ——_ + +...}{-+—- + ——_ +... 
+( ” 1! 82" sie +7 3! (8: a -t! a3 7 os! | 
(11) 
and fifthly the formula 
{ ee) \de = 2(1,—1,4+-1,—L4+4-—...). (12) 


0 

Miiller’s table was first of all extended to x = 20 by formula (8) and was 

later recomputed from x = 0 to x = 20 at intervals of x = 0-1 using 

formula (12), the results being given to 8 figures. This method depended 

on a table of I(x) which is as yet unpublished, and the author is indebted 

in this respect to Dr. J. C. P. Miller for kindly lending him the proofs of 

B.A. Mathematical Tables, Bessel Functions, Part Il. The results were 

checked by comparison with the first table as far as possible and also by 

differencing. The values of the integral for x = 20(1)25 were computed 

using formula (9) and as the differences are rather large in this region, 

zx 

a table of e-* I,(x)dx is given for 2 = 15(1)25. The results have been 
0 

checked by differencing. The differences 8?, 


= §2—0-184(84— 185); yt = 84/1000, 


and y*, (3) and (4), where 


8 


m 


are given for use with the formula: 


= E2 $9 2 2 4 
fo = (1-9) fot Of + £6 Sino + £5 8in + Mo 6+ Mi Hi- 

I wish to express my thanks to Professor R. O. Street for the interest 
he has shown and the helpful suggestions he has made in connexion with 
the computation, and to a referee for suggesting the use of formulae 
(10)—(12). 
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Interpolation formula: fs 





0 8668 6 24 3145 
= 3469 8150 26 6576 2 
2 3781 4628 29 2290 3 
3 4122 3859 32 O521 3 
4 4495 4119 35 1514 3 
55 4903 6452 6 38 5548 3 
6 535° 4946 42 2913 4 
wf 5839 7028 | 46 3949 4 
8 6375 3801 50 9017 4 
9 6962 0406 55 85135 
6-0 7604 6420 6 61 2880 
‘I 8308 6299 67 2603 
2 9079 9863 73 8209 
3 9925 2526 6 81 0282 
“4 1085 1738 5 8 8947 
5 1186 7285 5 9 7649 
6 1298 0639 10 7209 
7 1420 1376 Il 7715 
8 1554 0020 12 9265 
‘9 1700 8138 14 1951 
7:0 1861 8439 5 15 5901 
‘I 2038 4895 I7 1230 
2 2232 2861 | 18 8083 
3 2444 9217 | 20 6604 
4 2678 2517 22 6974 
75 2934 3162 5 24 9360 2 
"6 3215 3577 27 3975 2 
7 3523 8418 | 30 1041 3 
8 3862 4796 | 33 0802 3 
| 
‘9 $234 2521 | 36 3526 3 
8-0 4642 4372 5 | 39 9511 4 
‘I 5090 6395 43 90d5 4 
2 5582 8232 | 48 2611 4 
‘3 6123 3479 53 0470 5 
*4 6717 0083 58 3131 5 
8°5 5 64 1036 6 
6 7° 4738 = 6 
4 77 48047 
‘8 9736 7106 5 85 1890 8 
‘9 1068 6522 4 g 3668 
9:0 1173 0158 4 10 2998 
‘I 1287 6964 II 3260 
2 1413 7220 12 4554 
3 1552 2238 13 6977 
4 1704 4403 15 0651 
95 1871 7591 4 16 5695 
6 2055 6691 18 2245 
se 2257 8343 20 0465 
8 2480 0796 22 0510 
> 2724 4130 24 2573 
10°0 2993 04454 26 6854 
2 1 W252 4,4 4,,4 
mot EX Om +Miys T Mi yi 








z 


10” i} I(x) dx 


0 





au WN 


6 6 wo 


pwn 


6 SC wm AU 


— — - = -_ 
a — w we N 


OC COUVAU FHUKHHKHS CHOYAUW FHWNHHS CHIAW BWDH 


_ 
wm 


2993 
3288 
3613 
397° 
4363 
4794 
5269 
5792 
6366 
6998 


7994 
8458 
9300 
1022 


1124 ; 


12260 


£5 
1359 
1495 
1044 
1808 
1989 
2188 


2407 


4696 


=168 
5687 
6258 


6887 


1979 


2179 
2399 
2641 
2908 


2902 
3202 


3526 


0445 
4063 
I751 
2972 
o154 
8989 
8768 
2738 
8507 
8484 


5°54 
7159 
4183 
2760 
1210 
97°9 
0474 
7976 
g162 
3714 
4321 
6088 


373 
1156 
4453 
4261 
8961 


2861 
6802 
5882 
3489 
3973 
7795 
6363 
2429 
9916 
4196 
2206 
2607 
5960 
4919 
4440 


2045 


~ 


6854 
3575 
2989 
5362 
0g9g2 
0218 
339° 


5859 


9550 


46939 





I 389 
0980 


6506 
3413 
2042 
2568 
5182 
o1oo0 


7501 


7813 


z 
10" f I(x) dx 
0 





x» 


zs 


aint +> > ¢ 


2) 


4 


Ww Ww 


> + 


1652 
1819 
2004 





4475 8599 








1143 
9066 
5 0863 


0947 


9923 


6293 


O5oI 
8937 


4706 


4142 





— 





oe ¢ 








© | Lo(x) dx 





6742 

















ON THE HEAVING MOTION OF A CIRCULAR 
CYLINDER ON THE SURFACE OF A FLUID 


By F. URSELL 


(Department of Mathematics, The University, Manchester) 


[Received 20 April 1948] 


SUMMARY 

We consider the motion of a fluid of infinite depth which arises when a horizontal 
cylinder of circular cross-section oscillates with small amplitude about a mean 
position, in which the axis of the cylinder is assumed to lie in the mean surface. It 
is further assumed that the resulting motion is two-dimensional ; this assumption is 
justified when the cylinder is long compared with a wave-length, or when the fluid 
is contained between vertical walls at right angles to the axis of the cylinder. 
Expressions are obtained for the wave motion at a distance from the cylinder, and 
for the increase in the inertia of the cylinder due to the presence of the fluid. 


Introduction 

A CYLINDER of circular section is immersed in a fluid with its axis in the 
free surface. If the cylinder is given a forced simple harmonic motion of 
small amplitude about its initial position, a surface disturbance is set up 
in which waves travel away from the cylinder, and a stationary state is 
rapidly attained. When the cylinder is very long or when the fluid is 
contained between vertical walls at right angles to the axis of the cylinder, 
the velocity component parallel to the axis of the cylinder vanishes and 
the motion is two-dimensional. It is well known that at a distance of a 
few wave-lengths from the cylinder the motion on each side is described 
by a single regular wave-train travelling away from the cylinder, and 
that the wave-amplitude is proportional to the amplitude of oscillation 
of the cylinder, provided that the latter is sufficiently small compared 
with the radius of the cylinder, and that the wave-length is not much 
smaller than the diameter of the cylinder. 

In this paper it will be shown how the fluid motion can be calculated 
when the cylinder is oscillating vertically. The foregoing assumptions 
will be made, and viscosity and surface tension will be neglected. Then 
a velocity potential and a conjugate stream function exist, and it will be 
assumed that terms involving their squares may be neglected. From the 
potential or the stream function it is easy to deduce the wave-amplitude 
at a distance from the cylinder and the added mass of the cylinder due 
to the fluid motion. 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 2 (1949)] 
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Formulation of the problem 

Take the origin of rectangular Cartesian coordinates at the mean 
position of the axis of the cylinder. The z-axis is horizontal and per- 
pendicular to the axis of the cylinder, the y-axis is vertical, y increasing 
with depth. Define polar coordinates by the equations 

z= rand, y = reos8@. 

Since the motion is symmetrical about the y-axis, it is sufficient to con- 
sider the quadrant 0 < 6 < 47. The velocity potential ¢ satisfies 


24g 

: ta ar ¢ = 0, (A) 
ox* © oy? 

the stream function u satisfies 
ete nth (B) 
Ci” oy” 


On the free surface the pressure is constant, whence to the first order 
> od 1 . 
K¢- 0, 0 Aor, r>a, (C) 
oy 
where K = o?/g and 2z/o is the period (cf. ref. 1). Also, by symmetry, 
0d/00 0, 6 Q. (D) 
It remains to express the boundary condition on the cylinder. This is 
that the velocity component normal to the boundary just inside the fluid 
is equal to the corresponding component of the velocity of the cylinder. 
Suppose that the ordinate of the axis of the cylinder is 
y = leos(ot+e). 
Then at (asin «, acosa+lcos(ot+e)), the normal velocity is 


lL ods d y 


= — COS X, 
a Cx ( 
dy . 
ub — ©Y sin x, 
dt 


and to the first order this condition holds at (asin «, a cos a), whence 


— 


% = loasin(ot+e)sin@ on r=a. (E 
It is required to find a velocity potential and a stream function satisfying 
the boundary conditions (C), (D), and (E), and representing a diverging 
wave-train at infinity. To this end a series of non-orthogonal harmonic 
polynomials will be constructed satisfying (C) and (D); these will be 
superposed to satisfy (E) for values of Ka less than 37/2 by a numerical 
process. It is shown that when Ka is less than 1-5, this is permissible, and 
it will be supposed that the process is permitted in the wider range. 
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Construction of polynomial set 


It is easily verified that the set of stream functions 


2m 
a2" 





of (m = 1,2, 3....) 


| 


sin2mé, K_- sin(2m- m) cos 


y2m 2m—1 y2m-l sin 


is such that the conjugate velocity potentials satisfy (C) and (D), while 
on 7 = a it takes the values 


: ne... 
sin 2m@ + =—— sin(2m—1)@. 
2m— 


[t is clear on physical grounds that this set is not closed on r = a, since 
the sum of functions of the set tends to zero as r tends to infinity, whereas 
in fact the stream function for large r must represent a diverging wave- 
train. It is therefore necessary to add a function satisfying (C) and (D) 
and representing such a train of waves, e.g. the function describing a source 
at the origin (cf. ref. 3) 


gb [‘Y.( Kr; @)cos ot + '¥,( Kr; 8)sin ot], 


i 


TO 


where b is the amplitude at infinity and 
WY (Kr; 0) = me-Kreos8 sin( Kr sin 6), 


e—krsin@ 
Y( Kr; 0) = [ all sin( kr cos 6) 4 K cos(kr cos 6)! dk — 


e 


0 
—7e —Krcos@ cos( Kr sin 6). 


The functions of the closed set must be superposed so as to satisfy (E), 
and since there are no singularities on the cylinder, y% must be continuous 
on r = a when 0 < 6 < }r. 

Suppose, then, that the stream function % is expressed in the form 


mop ‘Y,( Kr; @)cos ot +-‘¥,( Kr; 8)sin ot +- 
gb . 





= 7 sin 2m@ K ~ sin(2m—1)0 
+ cos ot > Pom(Ka)a2™ , = = |+ 
: a 


2n—] y2m—l 


20 . a r . 
. “ - . o, |sin2mé K  sin(2m—1)0 
+ sin ot > Jo»(Ka)a2™ | —~ a. = ~- 
« y2zm 2m = r2m 1 
1 


» © r r 7 9 
where the coefficients p,,(Ka), qo,,(Ka) are assumed to be of order 1/m*. 
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This series converges uniformly outside and on r=a. On r=a the 
function ys is a multiple of sin 6, by condition (E), i.e. putting r = a, 
¥(Ka; @)cos ot + ‘¥.( Ka; 6)sin ot- 


8 


— 


1 


— oe Ka . 
| cos ot S Pom( Ka) sin 2m6 + a sin(2m—10| + 
, 2m 


+- sin ot S Jo»,(Ka) 


1 


C( Ka: t)sin @. 


zm 





sin 2m6 + a jsin(2m—1)0| 


To determine C(Ka;t), put 6 = 4m (say); then 
C(Ka;t) Y (Ka: 47r)cos ot 4 Y.(Ka; }7r)sin ot 
— , Re . 
+ cos ot > 15, (Ka)———- sin 4(2m— 1)z 
. Pam Sm l a ) 


Ka . 
sin $(2m—1)z, 





2m—1 


+-sinot S qo,(Ka) 
2 


whence it follows that z,,(Ka), q2,,(Ka) are the coefficients in the ex- 


pansions 
(Ka; 0)—¥,(Ka; 4x)sin 8 = > pon( Ka) fom(Ka; 4), 
1 
(Ka; 0)—¥,(Ka; 4z)sin 8 = ¥ qom(Ka) fom (Ka; 8), 
I 
where 
—_— Be ., , ‘ 
fon (Aa; @) sin 2m0 +-—— {sin(2m—1)@—sin @ sin 4(2m—1)r}}. 
2m : S 


Write 


i m-1lk 
Y(Ka;42)+ 5 \ rn Pan Ke) A(Ka), 
p hw 2m— ‘i 
1 
“ m-1k 
Y(Ka;4n)+ 5 2) Ae aan( Ka) = B(Ka). 
kt «3S 2M — = 
l 


Then the stream function on the cylinder is given by 
er Beebe i 
ys = | A(Ka)cos ot+ B(Ka)sin ot |sin 0, 
TO 
whence by comparison with condition (E) the ratio 


wave-amplitude 


aKa 


(42+ B)° 
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Calculation of the coefficients A(Ka) and B(Ka) 

The coefficients p»,,, Yom are the roots of an infinite number of equations 
in an infinite number of unknowns. For purposes of calculation, the system 
of equations was replaced by a system involving only a finite number of 
the polynomials f,,,,(Ka; 6). The functions 

Y.(Ka; 6)—'¥,( Ka; 47)sin 6, 

Y,( Ka; 6)—‘¥,( Ka; 32)sin 0 
were evaluated at = 0°(10°)90° by quadrature, with 

Ka/m = }, }, 4, 3, 2, 1,33 
and the polynomials f,(Ka; @),...,foy(Ka; 0) were fitted at these values by 
least squares, the corresponding coefficients being 


ian(Ka:N) (m = 1,...,N). 


The least squares condition provides a set of N simultaneous linear 
equations for p.,,(Ka;N) and similarly for q,,,(Ka;N). The matrix of the 
equations is symmetrical and the terms on the principal diagonal are 
larger than any other terms in the same row (except for the first row). 
The system could therefore be solved conveniently by relaxation methods. 
Trial calculations were carried out and it was found that an expansion 
in terms of six polynomials was adequate, giving a close fit at the chosen 
values of Ka and 6. Table 1 shows various functions of Ka to three signi- 
ficant figures. 
TABLE | 











mKa 

Kaln A(Ka) B(Ka) ,(A?+- B?) tan1(B/A) | a/(A?+ B?) m(Ka) 
° ° | 1°57 1°57 go | ° co 

1/6 75 | 2°23 2°83 —52 0°58 | 0-78 
1/4 2°75 — 2:06 3-44 37 o72 | 073 
1/2 5°64 +065 5°68 | +7 0°87 0°83 
2/3 6°17 4°39 7°57 35° 0°37 o°gI 

3/4 5°55 6°62 8°64 50° 0°86 0"94 

I 0°88 12°33 12°4 94 080 I‘oI 
5/4 12°6 II°2 16°9 138° 0°73 1°06 
3/2 —22°0 a7 22°0 183 0°67 1°09 








These calculations are based on six polynomials. The values o! 


,/(A?+ B?) and tan-(B/A) are given for convenience in interpolation, 
each of these functions being monotone increasing in the range; in fact the 
angle varies nearly linearly except near Ka = 0; when Ka = 0 
d . , 
———-tan-1(B/A) = 2 radians, 
d(Ka) 


as can be shown by expanding in a power series in Ka. 
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Virtual mass due to the fluid 

It is well known that when a long cylinder completely immersed in an 
ideal fluid is moving in any way perpendicular to its axis, the reaction 
of the fluid may be expressed in the form 


M’'f per unit length 


where M’ = zpa® is the mass of fluid displaced by unit length of the 
cylinder and f is the acceleration vector (ref. 1, Art. 68). This expression 
shows that the motion is unaltered if the fluid is supposed to be removed 
and a mass MM’ per unit length is added to the cylinder. The mass M’ is 
called the virtual mass due to the fluid. 

When the cylinder is in the free surface, the reaction of the fluid is no 
longer in phase with the acceleration. There is a component in quadrature 
which does a positive amount of work in each cycle and is thus simply 
related to the wave-amplitude. The component in phase with the accelera- 
tion does no work: it consists partly of a hydrostatic force which dominates 
the motion for small values of Ka, and which is easily seen to be 


2pa* d*y 


Ka dt?’ 
where y is the displacement of the cylinder. More interesting is the part 
due to the wave motion which will now be calculated. The velocity 
potential ¢ is easily derived from the stream function. It is given by 
70 


= * ®,(Kr; cos ot+®,( Kr; @)sin ot + 
qo 





2m—1 fs 


= , cos 2m K  cos(2m—1)0 
cos ot > Pom Kayan += = : " 


1 


, — _ . . feos 2mé K  cos(2m—1)0 
sinot S Jo»(Ka)a*™ | —,— +: - a 
Lat iis 2m—1 _ 
1 
from which the pressure —p ¢¢/ét can be derived. 
Here ©,(Kr;6), ®(Kr;@) are the harmonic conjugates of (Kr; 6@), 
(Kr: 6): 
®,(Kr; 6) me Kreosl egs( Kr sin 6), 


Z krsin@ 
emu fk cos(kr cos 6)— K sin(kr cos 6)} dk+- 
K2_/2 


0.(Kr; 60) 
+-ae-Kreos8 sin( Kr sin 6). 
It is seen that the pressure on the cylinder is of the form 


M cos ot+-N sin ot, 
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while the displacement of the cylinder is of the form 

P cosot+ Qsinaot. 
The component of the pressure in phase with the displacement is therefore 


MP+NQ 


’ PUG (P cos ot + Q sin ot). 


The force per unit length acting on the cylinder is 


40 
wa f 


. M, cos ot —N, sin ot), 


| pa cose dd (r =a), i.e. 
CL 


where 
47 


é __])\m-1 Gi 
M, [ ®,(Ka; @)cos 6 d@+ Zz Exe yr *Gam(Ka) inKaq,( Ka), 
: 1 


4m? 
0 
im . . - 
N, ® (Ka; @)cos 0 dé+ ba —. ~Pan(K) . 1 Kap,(Ka). 
: 4m*— 1] P 
0 


The displacement is 


l 
a (A sin ot— Beosot), 


aKa 
where A(Ka), B(Ka) are the functions shown in Table 1. 
It follows that the force component in phase with the acceleration is 


2pabg M, B+- NA 
a  A®1 B 


(A sin ot— B cosat), 


while the acceleration is 
bo? 


— —— (A sinot— Beosot). 
aKa 


The virtual mass is their ratio 


) 


- 


Mo B+NA 
A+ B? ° 


The values of the non-dimensional quantity 


M, B+N,A 
At, B 





m(Ka) = 


which may be described as an inertia coefficient, are given in the last 
column of Table 1. 
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It can be shown that as Ka + 0, 


m(Ka)—log —> 3—2log2—y —0-46, 

refore 1a 

vhere y is Euler’s constant. Hence the hydrodynamic inertia coefficient 
tends to infinity like log {1/(Aa)}, while the hydrostatic coefficient tends to 
infinity like z/(2Ka). 

The work done by the cylinder in a cycle must be equal to the energy 
transmitted by the waves at a distance from the cylinder in the same time. 
!), In terms of the various parameters, M, A—N, B = 3x*. This relation may 
be used to check the computations. 


Discussion of results 


The computations show that the amplitude ratio 


aKa 


,/(A?+ B?) 


is small for small values of Ka/z and increases as Ka/z increases, until 
Ka/z is approximately 0-6. As Ka/z increases beyond this value, the 
mplitude ratio decreases steadily throughout the range of computation. 
This effect may be ascribed to interference between waves originating 
from different parts of the cylinder surface. Qualitatively similar behaviour 
isexhibited by cylinders of various other sections. As Ka/z tends to zero, 


is 
the amplitude ratio tends to 2Ka. It may be shown that this result, which 
isin agreement with Holstein’s approximate theory (ref. 4), is still valid 
lor cylinders of section other than circular, provided that 2a denotes the 
horizontal diameter in the surface (beam). 

The inertia coefficient is of order 1 through the greater part of the 
range, except near Ka = 0, where it tends to infinity. On the other hand, 
the hydrostatic inertia coefficient tends to infinity even more rapidly, 
whence it may be concluded that in very slow heaving the deformation 
f the surface does not cause a significant change in the force on the 
cylinder. These results may be compared with the measurements made 
by Holstein (ref. 5) on a cylinder of rectangular section. Holstein measured 
the amplitude of the waves for various mean depths of immersion of the 
oweredge. When the depth is equal to about half the beam the amplitude 
should be comparable to that due to a circular cylinder. It is found that 

bial theory and experiment are consistent in the range of measurement 





0 < Ka/7 < 0-8). Holstein’s experiments on virtual mass (ref. 6) are 


too rough for comparison with theory. 
5092.6 ; 

| Q 

| 

| 
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Note on the magnitude of p,,,(Ka), do,(Ka) 
It has been shown that the calculation of the wave-motion requires the 
expansion of a function G(@), where 


0<0<}n, G(0)= G(4n) = 0, 


in a series of non-orthogonal polynomials 





. ° x ey ° ° 
fem(z; 9) = — sin 2m +- a {sin(2m— 1)@—sin @sin $(2m—1)z 


— 
a | 


me = 1.2.3.5) 
where x is a numerical parameter which is usually less than 5. 
Let the coefficients be denoted by 7,5,,(x) so that 
G0) = Ss Pom(X) fom(x; 4). 
m=1 
This expansion must converge uniformly throughout the range 0 < 6 < 4r; 
it has been assumed in the text that functions 7,,,(7) exist such that 
Pom(%) = O(1/m?) for fixed x. 
A proof that, in fact, 
Pom(x) = O(1/m?) (1) 
will now be given when 2 is assumed small (|x| < 1-5). 
THEOREM. Let G(6) be defined in the range 0 < 0 < 4r and let its second 


differential coefficient be of bounded variation. Then if |x| < 1-5 there exists 
an expansion 


G0) = & Pom (©) fom(5 8), (2) 

m=1 
5 “ , - A(x) . 
where | Pom(x)| <= (2m—1)3" (3) 


We first prove the following 


Lemma 1. (Limiting case x = 0.) 
Consider the expansion of G(0) in terms of the orthogonal polynomials 


Som(0; 8) = —sin 2mé. 


7 , - A(0) 
Then Pam)! < Se — 1s 
Proof. By Fourier’s theorem, 
47 
Pom() a. G(@)sin 2mé dé 


qi 


G"(0)sin 2mé dé, by integration by parts. 


7m 


9 
re 


O Oe ry SC — 
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Since G”(6) is of bounded variation, the last integral is of order 1/m (ref. 6 
§ 9.41), whence the lemma follows. 


> 
Suppose now that each coefficient p.,,(7) can be expanded in a power- 
series In 2 


De,(t) = > a@a* (m = 1,2, 3....). (4) 


n=0 


Substitute in equation (2) and equate coefficients of x”; then 


> aS? sin 2mé —G(8) (5) 
m=1 
J an» .. ‘ - 
¥ a” sin 2m0-4 S <™ ___ §sin(2m— 1)@—sin 6 sin 4(2m—1)z} = 0 
i a 2M—1 ni 


m=2 


(n>1). (6) 


Suppose further that the infinite series in equations (5) and (6) converge 
uniformly throughout the range. It is then permissible to multiply each 





equation by sin 27@ and to integrate term by term. 
From equation (5), 


$or 
Lrayy. | G()sin 276 dé. (7) 
0 
From equation (6), 
Dn b (n—1) 16 2 
Leva”) = Y Bm (—1)m+r _(2m— ae (8) 
‘ 4r°@— 1 Ly 2m—1 (2m—1)?—(2r)? 


m= 


It must now be shown that equation (8) defines a double array a; that 
the corresponding power-series defined by equation (4) are convergent for 
small x; that these power-series satisfy (3), and that 





Y fom(vs0) Y agra" = G6). 
m l n 
LEMMA 2. 7'o show that 
A(0 
ag) < <0) | (9) 
Pp (2r—1)8 


Proof. From equation (7) a{?) = ps,(0), so that the result follows from 
Lemma 1. It also follows that equation (5) holds. 


LemMa 3. 7'o show that 


A 
* ; where A 


P<. n = (§)"A(0). (10) 
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Proof by induction on n 


From Lemma 2, equation (10) holds when n is zero. Suppose that 


equation (10) holds for n = 0, 1, 2,..., N—1. From equation (8), 
x N-1) 2 
jray)| = i| 2, a Lyne (2m — 1)? é 
2m- Paes, 
or . aN - 1) (2 m- aa 2 
4r2—1 a, (2m—1)|(2m—1)2—(2r)? 
212 Ay 1 S ] _ 
4r2—] —, (2m— 1)?| (2m—1)?—(2r)? 
_ WAyial l 
~— 4r2—1 | Ly (2m—1)*{(2r)2?—(2m—1)*3 | 
412 —_ - 
) har (2M—1)?{(2m—1)?—(2r)?! 
Now 
A ] 
>. (2m— 1)?f(27r)2— (2m— 1)?! 
lfx< = l 
4r?| La (2m—1)? ‘3 (2r)?—(2m—1)? 
4 m 
l ] 
al eo) 
ipl ) L6r* 
Che second series may be written 
. | 
a ; (2m — 1)?{(2m— 1)?—(2r)?! (2r-4 1)2(4r +1) © (2r+3)*(12r- 


l = 
© (2r+5)?(20r+25) 5 | p3 , (2m—1)?f(2m—1)? 
asf 


Now 7 ] 


)2(127r+-9) ~ 400 


yr? ] 
2 


(2r+5)2(20r+5) ~ 1000 


— 
bo 
ww 





9) 


(2r)?} 
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< l [ du 
L,  (2m—1)*%{(2m—1)?—(2r)*} J (2u—1)3{(2u—1)?—(2r)?} 
oes r+3 


| du 


2(20r+- 25) | (2u—1)8 
r+3 


l 1 
4(20r+-25)(2r+5) ~ 160r?’ 


< if trea tiel <i 
7 





- — — — . 
LZ, (2m—1)%4(2m—1)*—(2r)*! 45 400 1000 160 * 3lr? 
1 
Hence 
1a \ 2rA 7-3 d aie 2 
sill (4r2—1)r2\16 | 31 
als gAy _ Ay 
(2r—1)® (2r—1)8 
[his proves Lemma 3. 
It follows immediately that 
. Ady, J 4 bs P 
= 7 isin(2m 1)9—sin @ sin 4(2m—1)z' (11) 
Lag 2M - 


s uniformly convergent. Its associated Fourier sine series is 
> ay) sin 2mé (n 1), 
here aS” is defined by (8). It follows from Lemma 3 that the last series 
mverges throughout the range, and so its sum is equal to (11); that is, 
equation (6) holds. 


It also follows from Lemma 3 that 


Pom\") 3 Am x" 
n=0 
s defined for |a 1-5: again 
A(0) : — ; 
Do,,(2) 1-+-4|x|+(s)*|2|° 
2m ns 
A(x) ”" 
= (a 1-5) (12) 
Zn 1) 
ind since | f,..(a; 0) 1+2\a > Pom (©) fom("3 9) converges. 
m 1 


It only remains to be shown that 


> Pom(X)fom(2; 9) = G8). 


l 
m l 
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Lemma 4. Jo show that 





lim 4, £ Som: 13 6)( > adm x) = : > Pom(2 ) fom ;). (13) 
N-o m=1 m= 
. | A(x) cot 
For | f. : aig an| <- ——. (1-+- 2/21); 
Som (%: > 2m | (2Qm— 1! | ) 


e = A(ax)(1+2|a]) - 
since 2 ps COnverges, the series converges uniformly with 
2m— 1) , : 

m=1 ( 
respect to NV, by Weierstrass’s M-test (ref. 2, § 3.34); hence equation 


(13) holds. 


Proof of the theorem. Consider the expression 


Lo 8) 
(n) 
y Som > Agm & 
m=1 n=0 
From (5) and (6) 
0 
; ae (n) y 
> Fam( 2s 8)( > Agm v n)— -G(@) 
m=1 
ws a? 
—yN+1 3 —2™ _ fsin(2m—1)0—sin }(2m—1)z sin 6}, 
9 - J 
; 2m—1 
m=2 
but 
ln pN+1 s a{sin(2m— 1)@—sin 3(2m—1)z sin 8} 
m=2 
< 2|a|N+1(2)V.A(0) > — (from Lemma 3), 
(: 2m— 
which tends to zero, as NV tends to infinity provided that |2| < 1-5; so 
that . 
lim s Samle 6) > aa G(@). 
N-o m= = n=0 


x 


From Lemma 4, the left-hand side is ¥ py,,(7)fo,,(2; 9); also from (12) 
m= 


; A(z) 
Dante)! < aud 
Pam\*?! (2m—1)8’ 
which proves the theorem. 
Note. The foregoing proof applies only when |a| < 1-5, but it has been 
assumed that the theorem is valid in a larger range including |x| < 1°5. 
[ am indebted to the Mathematics Division of the National Physical 


Laboratory for the computations in this paper, and to the Admiralty for 
permission to publish. 
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THE BALLISTIC EFFECTS OF BORE RESISTANCE 


By J. CORNER (Armament Research Establishment, 
Fort Halstead, Sevenoaks) 


[Received 7 July 1948] 


SUMMARY 

Resistance to motion of a projectile down the bore of a gun alters the pressure and 
velocity attained, but is usually neglected in ballistic theory. It is shown that 
weighting factors can be calculated for each point in the shot travel, such that the 
effect of a bore resistance with any desired dependence on travel can be quickly 
found by integration of the product of weighting factor and resistance. The method 
fails if the resistance is sufficient to stop the shot for a time, and initial engraving 
resistance is therefore not a special case of this solution. 


1. Introduction 
THE classical central problem of internal ballistics is the calculation of the 
pressure, travel, and time relations in a gun, given the characteristics of 
the charge. The number of solutions advanced is already very large, and 
when assembled, as in the compendium of Cranz (3), they tend to produce 
the feeling that little remains to be done. This is not correct. The initial 
resistance to motion,’ due to the engraving resistance of the driving-band 
of the projectile, is commonly taken into account, often as a ‘shot-start 
pressure’. On the other hand, it usually happens that resistance which 
occurs while the shot is moving down the bore is included, if at all, by 
assuming that this resistance depends on time in the same way as the 
accelerating pressure, so that the effect is exactly equivalent to an increase 
of shot mass. This is a very special form of resistance, a poor approxima- 
tion to the types found in practice,} and while it does not enable the shape 
of the resistance-travel curve to be varied in studying a given firing, yet 
the form changes against one’s wishes when the charge is altered. Among 
published work the only solutions of greater value are those of Proudman 
(4), Voituriez (6), and Clemmow (2). The two former writers assumed 4 
constant resistance, which is a special form but has the virtue that it does 
not change its shape with the ballistic parameters in the solution, and is 
often approached in practice. These authors found the essential parameters 
which specify the solution, though they did not give any explicit results, 
which would apparently have to be found by numerical integration. 
Clemmow studied the effect of bore resistance with the approximations 
of a covolume equal to the volume of the solid propellant, and y, the 
ratio of specific heats of the propelling gas, equal to unity until the charge 
+ Except that it is often a fairly good approximation in small-arms firing jacketed shot. 


(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 2 (1949)] 
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Clemmow found 
the most general law of resistance which allowed a solution in finite terms; 


is completely burnt (the ‘isothermal approximation’). 
he studied also the case of a constant resistance, finding a form for the 
solution which reduced the labour of computation of any particular 
example, and an approximation for small resistance. The Clemmow 
resistance law covers a fairly wide range of types of behaviour and shows 
in an elegant way the nature of the effects produced by such kinds of 
resistance. This solution is, however, not easily applied in practice when 
one is given the details of the gun and an actual resistance-travel curve; 
the determination of a Clemmow curve is laborious and it is not always 
possible to come reasonably close to the assigned resistance law. There is, 
of course, also the restriction that y = 1, which has been superseded in 
most modern ballistic systems. 

In the unpublished literature ballistic results taking bore resistance into 
sccount have been computed by numerical integration for a number of 
special forms of resistance. For resistance depending on travel in any 
other way it would be necessary to repeat the entire computation. 

Thus there has been no simple rapid way to find the effect of a bore 
resistance depending in any desired manner on the travel of the projectile. 
Such a method is given in the present paper. As an introduction it may 
be helpful to state, in a general way, the experimental results which led 
to this work and some of the technical problems to which it can be applied. 

Evidence for a resistance persisting for a substantial part of the travel 
of the shot has accumulated from many directions; much of this evidence 
is only qualitative, though its bulk and its general consistency are impres- 
sive. It is found that the fitting of theoretical ballistics to observed results 
often leaves a discrepancy which could be attributed to a mean resistance 
equivalent to a pressure of order 1 ton/sq. in. It is true that in many cases 
one might reduce the discrepancy by altering the details of the theory or 
by postulating small errors in some of the data assumed; however, such 
ad hoc manipulations are, in bulk, less plausible than the idea of a bore 
resistance. From the experimental side, the strains recorded in the barrel 
ive sufficient to indicate an appreciable axial resistance even when the 
coefficient of friction has fallen as low as 0-01. A resistance is also usually 
indicated when attempts are made to correlate measured space-time curves 
with pressure-time curves recorded at gauges spaced along the barrel. 
Considerable progress was made in this field during the war, though the 
only published work is that of Rossmann (5). Pre-war work was sum- 
marized by Bodlien (1), and was mostly on small-arms. Here the whole 
body of the projectile is engraved, being kept pressed against the rifling 
by the deformation under the high acceleration in the gun; a substantial 
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resistance to motion is therefore to be expected in small-arms. Unfortu- 
nately the time-scale of the phenomena is shortest in these guns, and it 
has not been until recent years that the natural frequency of the recording 
system has been raised sufficiently to give significant results for the 
resistance. 

It can be taken, then, that in many cases a resistance does persist over 
much of the travel of a projectile, and it is desirable to investigate its 
ballistic consequences. It is clear, too, that knowledge of the resistance is 
not yet in a final form, and that, therefore, it is better not to restrict 
attention to any particular form for the dependence on travel. For many 
purposes it is sufficient to be able to ‘touch-up’ some ballistic solution 
(which neglects bore resistance) by adding a correction to be computed 
for the particular resistance in use. 

In this paper a method is given by which the changes in muzzle velocity, 
peak pressure, and position of ‘burnt’ can be computed quickly for a bore 
resistance of any desired form, provided only that the resistance is nowhere 
sufficient to halt the shot for a time. The change in the complete pressure- 
space curve, which is of interest in gun design, could be found from the 
same formulae, but it is quicker to make an accurate numerical integration 
after the desired peak pressure or muzzle velocity has been obtained by 
using the results of this paper. The condition that the velocity of the shot 
must not fall to zero during the application prevents initial resistance to 
motion being included as a special case of our bore resistance. 

The line of approach is the choice of a standard form of resistance, the 
determination of its ballistic effects, and the construction of more realistic 
resistance laws by superposition of standard forms suitably spaced along 
the travel. Analogous methods are common in applied mathematics. In 
external ballistics, for example, the change in range due to perturbations 
such as head wind, which occur with varying strength at different heights, 
are obtained by integration over the actual distribution of wind, weighted 
by a factor which represents the effect of unit wind concentrated along 
unit length of the trajectory. This ‘weighting factor’ for the effect of head 
wind on range is a function of position along the trajectory, and is one 
of a series of analogous weighting factors which are computed from 4 set 
of equations adjoint to the primary ballistic equations. 

A similar practice is found in operational methods for sets of linear 
differential equations. The behaviour of the system is investigated for a 
standard impressed force, the Heaviside unit function. The response to 
a force varying in any way with time is obtained by splitting this force 
into a sum or integral of unit functions suitably spaced in time. The 
response to the complex force is, by the linearity of the system, equal to 
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the sum of the responses to the individual unit functions. Since the 
equations of internal ballistics are not linear, it cannot be taken for granted 
that the ballistic effects of a set of ‘unit bore resistances’ can be super- 
posed ; that the error is small has to be verified. There is one similarity 
between the Heaviside unit function and the unit bore resistance which 
will be chosen in this paper: the usefulness of the former is in no way 
impaired by its artificial character (for it is an idealization which can be 
approached but not attained in practice); likewise it is possible to use a 
standard bore resistance which is an idealized form, so long as by super- 
position one can construct any bore resistance likely in practice. 

The method given in this paper has been applied to some technical 
problems and abnormal behaviour of guns, thought to be due to bore 
resistance. One may mention certain effects produced by chromium 
plating of machine-gun barrels, or by using tapered bores, the effect of 
grease in a gun on the ballistics of the first round of a series, and the 
variations of muzzle velocity during repeated firing. 

It is to be expected that when the bore resistance is better understood 
it will be found desirable to choose simple forms and to compute ballistic 
tables with such resistances taken into account in an accurate way. This 
would be a task of some magnitude, particularly if a wide variety of shapes 
had been found to be necessary, and though the results of the present 
paper form only an interim stage, it is possible that for a long time this 
may be all that is justified by our knowledge of bore resistance. 


2. Equations of unresisted motion 

Let A be the area of cross-section of the bore. The standard bore 
resistance will be taken as a constant, AP,, over a short travel Az near the 
point 2), and zero at all other shot positions. The effect on the ballistics 
can be calculated, to the first order in P,Ax, by subtracting an energy 
AP, Ax from the kinetic energy of the projectile as it passes through the 
point x». The ballistic equations from the point 2) onwards are exactly 
the same as before except for a perturbing term linear in P, Az. If the 
ballistic system is a simple one this perturbation can be integrated exactly. 
If not, numerical integration would be necessary, in which case it would 
probably be quicker to abandon weighting factors altogether. 

One is thus led to consider the choice of a suitable ballistic system to 
which to apply this method. As has just been seen, it must give an 
analytical solution, the simpler the better. A second point can be revealed 
by letting the resistance of length Ax occur at shot-start. Then, if P, is 
finite, there is a finite change of muzzle velocity even if Ax tends to zero. 
The weighting factor which multiplies P, Ax must therefore tend to infinity 
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a8 % moves towards shot-start. Hence the weighting-factor method fails 
to represent initial resistance at all, and must become increasingly 
inaccurate as the position of the resistance approaches shot-start. It js 
therefore immaterial whether the basic system, before the addition of 
bore-resistance terms, be a good or bad representation of initial resistance, 

In this paper there is used a simple ballistic system, well known in the 
open literature. It is assumed that 7, the covolume of the propellant gas, 
and 6, the density of the solid propellant, are related by » = 1/8. This 
‘neglect of covolume terms’ simplifies the equations considerably, but at 
the same time it is possible to make partial correction for the error by 
changing the ‘force’ of the propellant. 

Conventional assumptions are made about the Lagrange pressure- 
gradient in the propellant gas. 

The burning law assumed is 


@ = (l—f)(1+ Af), (1) 
df 
J = = Pr. 2 
f dt - \ 


where ¢ is the fraction of the charge burnt up to time ¢, D is a length 
characteristic of the effective size of the propellant grains, 8 and @ are 
parameters, constant during any one firing, and P is the pressure at the 
breech. The variable f could be eliminated from (1) and (2), and, indeed, 
this is often done in continental treatments of internal ballistics; f is, 
however, a useful auxiliary variable in the process of solution. From (1), 
it follows that f 1 at the start of the firing, where ¢ 0. 

An ‘unsophisticated interpretation of (1) and (2) is that 4fP is the rate 
of recession of a cordite surface at pressure P, and that the grains are long 
compared with their cross-section, which is a simple shape with fD as the 
least distance between exposed surfaces; @ is then obtainable from the 
geometry of the grains. For example, long cords would have @ = 1, Das 
the initial diameter, and f/D as the remaining diameter at time ¢. Long 
tubes would have @ = 0, D as the initial annulus, and fD as the remaining 
annulus. In practice there are many factors which alter the interpretation 
of (1) and (2) while not removing their empirical value. The rate of burning 
of a cordite surface often increases less slowly than P and is also influenced 
by transverse gas velocities, which may be important in shapes with 
perforations. Initial resistance to motion is often counterfeited by increas- 
ing 8. The heterogeneity of any practical charge increases 6, as also does 
non-uniformity of any kind within the single grain.+ Short chopped grains 

+ For the effect of an eccentric perforation in tubular propellant, cf. Kleider, Zeit. ge. 
Schiess- und Sprengstoffwesen, 38 (1943), 131. 
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or cubical or spherical grains do not satisfy (1) exactly, even on the 
seometrical interpretation, but these propellant shapes are of little impor- 
tance. It remains true that (1) and (2) lead to simple formulae and most 
of the difficulties can be circumvented, at least to practical accuracy, by 
suitable choice of 8 and @. 

Let V(t) be the velocity of the shot, 2 its travel, U the chamber capacity, 
Athe bore area, C and W the weights of charge and shot, RT, the ‘force’ 
of the propellant, y its ratio of specific heats, suitably increased to allow 
for heat loss, and let 7, 7’(t) be the gas temperature, averaged throughout 
the volume behind the shot at time ¢. A correction for the rotational 
inertia of the projectile can, if desired, be included in W. The temperature 
of the propellant gases when formed without performance of external work 
is 7), and R has the meaning conventional in ballistics, namely, (gas con- 
stant per mole)/(mean molecular weight of the propellant gases). 


The equations of unresisted motion are 


d (] J)G ; Of), (1) 

df 
rE pe or’. 2 
dt (*) 

dV 
(Vi 1C Af. 3 
. hi s) 
1 P(x+l) CRT, 4T'(1+4C W)/- 4C/W), (4) 
4(1—T") = (y—1)(W+40)V2/2CRT,, (5) 

where the length / is defined by 

Al Uf —= C78: (6) 


These equations are orthodox and do not need any commentary. 





3. Solution with standard bore-resistance 
For the standard bore-resistance one may take a thrust AP, acting 
through a small distance Ax at a point x), where f = f,, V = Yj, and so 


| on.t The factor A is included for ease of visualization of the resistance 
| as an effective pressure. It is supposed, to begin with, that Az is infini- 


tesimal and P, Ax finite. 


This bore resistance acts for a time Aa/)), which is assumed to be 
infinitesimal. The only immediate result is therefore a decrease of the 
| kinetic energy of the shot by AP, Ax. Take first the case where the 
| resistance occurs before ‘burnt’. For travel less than 2, the solution is 
the normal solution of equations (1)-(5). For travel greater than 2, only 
the energy equation (5) needs alteration, though there has also been an 
| instantaneous drop in the velocity V. 


No confusion need be feared with the ‘force’ RT,. 
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One should note two limitations of this method. Firstly, if Y% is zero 
Ax/V is not infinitesimal; hence initial resistance cannot be represented by 
this model. Secondly, if the kinetic energy of the shot wken it meets the 
resistance is less than AP,Ax the method breaks down, for it gives the 
projectile an imaginary velocity after crossing the resistance. In practice 
the projectile would come to rest and would restart after the pressure had 
built up behind the shot, that is, after a finite time. This is again a matter 
of Ax/V, not being infinitesimal. 

Just before meeting the resistance the kinetic energy of the charge and 
projectile is 4(W+4C)V?; immediately afterwards the kinetic energy is 
3(W+4C)V?—AP, Ax. The instantaneous drop of velocity of the pro- 
jectile is transmitted through the gas with a finite velocity, namely a,, 
the velocity of sound in the propellant gases, and the time taken for the 
waves to settle down is of order 2(U + Ax,)/Ady, which is comparable with 
the time taken for the shot to reach the muzzle. The problem of the 
pressure-distribution in the column of propellant gas is therefore more 
difficult than in a conventional solution. If the velocity of the projectile 
immediately after the resistance is written as V,, then it is assumed here 
that the kinetic energy of the system at that moment is 4(W+4C)V?. 
In reality it should be 4WV?+-34Cv?, where v(t) drops from Vj, towards the 
V(t) curve. The approximation here adopted replaces this gradual 
approach by an immediate transition of v from ) to V. The extent of the 
error which could arise is small except for guns with C of order W, and 
velocities therefore around 4,000 ft./sec., but then the conventional 
Lagrange corrections also fail for such guns. 

With this approximation, then 


(W+4C)V? = (W+4C)V2—2AP, Az. (7) 
Equations (1)—(4) remain true after the resistance, and (5) is replaced by 
2C RT) d(1—T") = (y—1)[(W+-40)V?2+2AP, Az]. (8) 
From (3), V =V,+AD(fy—f)/B(W +40), 0 
so that (8) becomes 
20RT, $(1—7") = (y—1)(W+40)[V3-+{AD(fo—f)/B(W+40)}?+ 


+2ADVK(fo—f)/B(W+40)}- (10) 
Combining (4) and (10), 


2CORT,$ = (y—1)(W+4C)V3+{AD( fo—f)/B(W+40)}?+ 
+24 DVi(fo—f)/B(W+3C)]+24 Pew +-l(W+40)/(W440), (1) 
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which may be written as 
1—f )?— 2 ae Y)(1 —ho) Mfo—P) 
2B W +40) 
-AP(@-+1)(W+40)|(W+40). 








ORT, 4 = (y—1)(W+40C) 42D"! 


Up to this point the equations have been exact. Now, keeping only first- 
order terms in (7), (9), and (11): 


V, = %—AP, Aa|(W+40)h, (7a) 
V = AD(1—f)/B(W+4C)—AP, Aa/(W+4C)h, (9a) 
20RT)S’ = (y—1)(W+-4C)[{AD(1—f)/B(W+4C)P— 
2AP, Axl fo—f )/(1—fo)(W +40)]+24 P(w+l)(W+4C)/(W +40). 
( 


lla) 
It is convenient to introduce the ‘central ballistic parameter’ 


M A2D?(W +-4C)/B°?C RT (W+4C)? (12) 
so that (lla) can be written as 
(1—f)(1+6f) = (y—1)M(1—f)?/2—(y—N) AF, Ax(fo—f )/(1—fo) C RT) + 


(a4 FLAP dal (1—f,)CRT,—M(1—f)], 
Cc 


leading to 


Z4+-M(a ne AP, Aa (x Fife] /—Na—-fJORNy 

(13) 
where Z = 1—}(y—1)M+4+6,f (14) 
and 6, = 0+Hy—1)M. (15) 


Equation (13) shows the first-order effects of the resistance on the equation 
of unresisted motion 
, df 
Z+-M(x+l) = 0 


16 
= (16) 


and one can obtain a solution correct to terms of order P,Ax by using 
(16) on the right of (13). This gives 


2+M (e+) -AP, Aal Z/M-+(y—1(fo—f A—S )(1—fo) CRT), 


and so | 
dex tdf, AP,Ax fo—f) af 
- i) a ob 0 _ —1)] mE 
| wap’ rf Z° CRT,(1— | | ao Zi Lot ad oreratt 
(17) 
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The lower limit in these integrals corresponds to the point, suffix zero, at 
which the resistance occurred. The integrals are easily evaluated: 

0, { df/Z=Im(Z/Z,), 9, | df/Z2 = 1/Z,—1/Z, 

0 0 
) | df/Z.i—f) = In{Z(1—fy)/Z,(1—f )} 
0 
(1-6)? [ df/Z21—f) = In{Z(1—fy)/Zo(1—S)} + (1 +-0)1/Z—1/Z}. 
0 


Substitution in (17) gives 


In(1-ta/l) = “ing +.0)/2}4+ 


6, 
*- AP, Ax {] +-@— (y—1)J ~~ to) In (Z(1 fi So) |. 
CRT(1—fy) ye "Zell =p 


+ (y—1)M{1+0—80,(1—f,)}{1/Zp—1/Z}/0,(1 +8)]. (18) 


Solution with ‘unburnt’ at ‘muzzle’. In this case the value of f at the 
muzzle is computed from (18) and the muzzle velocity follows by (7) and 
(9). This case is not of practical importance and one need not test the 
accuracy of superposition and higher-order terms, which may be expected 
to be of the same order as when ‘all-burnt’ occurs inside the gun. 

Solution for ‘all-burnt’ inside gun. This is the important problem. It 
will be remembered that one is dealing here with a resistance during the 
burning period. Resistance during the adiabatic expansion will be dis- 
cussed later. 

One needs the change in In(x,+-/), where suffix B refers to conditions 
at ‘all-burnt’. This change is 

hints) =. E +O (y= VM (So) 17,1 —f,)| Zo} 

CRT,(1—fy) (1+)? ay 
1) M{1+0—0,(1—f,)}{1/Z,—1/Z,}/0,(1 6)}. (19) 


Let suffix Z refer to values at shot-ejection. Then integration along the 
adiabatic from ‘all-burnt’ to the ejection of the shot gives 


Vi, = V3z,+AP,(2p+)0/(W+40C), (20) 
where ® = 2{1—r!-}/(y—1) )] 
and r = (t,+1)/(xp+l). (22) 


In (20), Vz, Pp(v@p+/), and ® are all altered by the resistance. 
From (11a), 
AP;,(%p+)/(W+4C) = CRT —h(y—1)M}/(W+-40C)- 


rly 1)AP, Aafy (1—fy)(W + $C). (23 
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Also A® = 2r-vAr = —2r!-YAln(x,+), (24) 
and from (9@) A(Vz,) 2A P, Ax/(1—fy)(W+40), (25) 
so that 
(Vz) 2A FP, Aa{1—fo+ for %}/(1—fo)(W+40)— 


21H MMAR, Ae (1+ 67 DMO Sane gf) + 


(W+40)(1—fy) (1+) 
+ (y 1) M{1+0—8,(1—fo)}(1/Z>—1/Zp)/0,(1+8)), 
and finally 
(W+40Vx(1—fy) AV ,/AP, Ax 
(1+0—(y—1)M(1 
(1+6)2 
| (y—1)M{1+0—6,(1—f,)}(1/Zy—1/Zz)/6 (1-+8)]. (26) 


This gives the change of muzzle velocity AV,, due to a bore-resistance 
AP, Ax occurring at the point (fy, Z)), in terms of the characteristics of the 


] fy for v¥—y" "Zp Soin {Z; 3(1—fy)/Z Lo} Ty 





solution without resistance. 
For resistance very early in the travel, f, is nearly unity, and 

l-yF | 
te nf +6) ) 

8)(1—fy) \Z "= (1—fy) )’ 
showing that for resistance ae adean early in the travel the 
muzzle velocity rises. For resistance at ‘burnt’, fy = 0 and 

q 1/‘'\/7 y > ae 
(W+4C)V,, AV,,/AP, Ax = —1, 

and the muzzle velocity falls. When resistance occurs during the adiabatic 
expansion the kinetic energy of the shot and gases at the muzzle are 


(W+40)V,, AV,,/AP, A ~ 


reduced by the work done on the resistance, without other effect on the 
ballistics (except on the time to traverse the bore). Hence 
(W+4C)V, AV,/AP, Az = —1 


if the resistance occurs after ‘all-burnt’. 


4. Accuracy of first-order theory 

The results of the preceding section were obtained by expansion in 
powers of AP, Ax, keeping only the first term. To find the importance of 
higher-order terms it is easiest to take a specific example. A heavy A.A. 
gun was considered, and a charge of 8 lb. 5 oz. of cordite, giving 20 
tons/sq. in. peak pressure, was assumed. With 6 = 0-3, y = 1-34, and no 
covolume term, the velocity from the basic method is 2,642 ft./sec. Peak 
pressure occurs at a travel of about 4 calibres and ‘burnt’ about 10 calibres 


farther on. A series of exact solutions, with exact allowance for the 


5092.6 
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assumed resistance, was computed on a National accounting machine, 
To get sufficient accuracy in the differences AV,, these solutions were 
computed to an accuracy of 0-1 ft./sec. 

Table 1 shows the effect of increasing P, Ax at a given point early in 
the travel. At this point, where 2, = 0-58 calibres, the shot would be 
brought to rest by an instantaneous P, Ax of 7-99 ton/in., so that one 
would not expect this theory to be more than extremely rough at such 
values of Ax. The agreement shown in Table 1 is therefore very satisfy- 
ing. The deviation of AV,,/P, Ax from first-order theory is only 7 per cent. 
when F, Az reaches half the value needed to stop the shot. 


TABLE | 





(W+3C)(1—fo)Vz AV eg 


P, Ax/7-99 AP, Ax Method 
° 0624 First-order theory 
0°25 0-654 Exact integration 
o'5 0°67 a 
I 0°85 ” 





Accuracy of first-order theory at 0-58 calibres travel (f, = 0-9). 
P, Ax = 7:99 ton/in. stops the shot. 

Another way to show the agreement is to calculate the effect of 
P, Ax = 7-99 ton/in. at various positions along the bore. Table 2 shows 
the comparison to be extremely satisfactory except when the resistance 
stops the shot. 


TABLE 2 





(W+4C)\(1—f.)VpAVg/AP Ax | 
fo | Exact integration | First-order theory | Error 





o’9 085 0624 | 0°23 
od 0°25 0°231 0°02 
o-7 — 002 0°026 | ool 
06 0°22 0°226 | O°or 
04 0°53 0°540 | oor 
o°2 —o78 -0:789 | oror 
fr) —I —I | 0 
Standard bore resistance, P, Ax = 7-99 ton/in. at various positions along the bore. 


One may conclude that first-order theory is adequate to represent the 
effect of a concentrated bore resistance except when this is sufficient to 
halt the projectile. 


5. Calculation of the effect of long stretches of bore resistance 
One often wishes to find the effect of bore resistance extending over 4 
large part of the travel. Such a resistance, no matter how it varies, cal 
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be represented as an integral over suitably chosen resistances of our 
standard type. To see the errors likely in this process it is sufficient to 
consider the special case of constant resistance, and since the exact 
solution is known from ‘all-burnt’ to ‘muzzle’, one need consider only the 
part up to ‘all-burnt’. The weighting-factor theory makes 


"AV, dx », 
AV; | ea df 
| AV; M(x ) af 
Ax Z . 
MA * («+P (1—f Ve AV;,(W +40) 7 
TW oe ef. (2 
Vie(W +30) | (J al AP, AX q (34) 


This integral is easily evaluated. P, is a known function of x, which is 
itself a simple function of f. So is Z, and the ‘weighting-factor’ in brackets 
comes from (26). Thus (27) is easily computed by a numerical integration, 
whatever the variation of the bore resistance. 

The examples solved exactly on a National accounting machine had a 


constant resistance equivalent to P, = 1 ton/sq. in. of bore area, extending 
from ‘burnt’ to f, = 0-5, 0-7, and 0-9. Table 3 shows that the agreement 


with (27) is very satisfactory. 


TABLE 3 





* Oseon t8 | AV (ft.[sec.) 





a | r : 
n in bore solution Exact | (27) Error 
to 13°43 calibres og too 23°3 23°5 | o2 
to 12°4 o'7 too 29°38 | —30°7 o'9 
» to 12°42 + os too 20°90 —27°5 | 06 
a 
The effect of a constant bore resistance for a considerable part of the travel. 


P. 1 ton sq. in. 
Another simple test is the effect of P, = 1 ton/sq. in. from 0-58 calibres 
travel to the muzzle. This is —115-9 ft./sec. by exact calculation, and 
116-1 ft./seec. by using (27) up to ‘burnt’ and thereafter the exact 





solution 










(W-+40)A(V3,) 2AP(%_p—2Xp). (28) 


6. The effect of bore resistance on peak pressure 
This can be treated by a weighting factor in the same way as muzzle 
velocity. The peak pressure is, however, much the less important of the 
| two, being needed only if one is ‘touching-up’ a ballistic system which 
| omits bore resistance. Furthermore, the peak pressure is affected only 


by resistance in the first few calibres of travel. 
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From (11a) and (18), 
In P = Inf[CRT,(W-+4C)(1—f)/(U—C/8)(W+ 10)]—F nc +0) + 
+ (1+ M/6,)In Z+- 





Aj Ax [ -1; (fof), 140—(—-M fo) j,{ 20-—P)) _ 
" ORT(1 —fo) he (1. + \Z( Ale 
—(y—1)M{1+6—6,(1—f)}(1/Zy—1 |Z)/0,(1+8). (29 
The maximum pressure P,, occurs when 
ee 
df 


and the solution is 


yM+6— 
f= yM +20 


The change in InP. 


ae term proportional to AP, Ax/C RT,(1—f,). 


»n 18 Obtained, to the first order in AP, Ax/C RT,(1—f, 
by writing f = (yM+6—1)/(yM-+-26@) throughout (29). Hence 


CRT,(1—f,)AP,,|P,, AP, Ax 


= M+ 20K aa +1—0—yM}/(1+-6){8+-4(y+ 1M} 
ae —fo)} : 
a ] 0) { 1 i 188 
° ( oF In[ Zo/(1—foX9+3(y + YM] 


+ (y—1)M{1+6—6,(1—f,)} > 
. . -8){8-+-4(y-+1)M}Zp]/0,(1-+0)(0-+-4(y+-1)M}. (30 
This holds if f, > (yM+0—1)/(yM-+26) > 0. If 
* < (yM+6—1)/(yM-+26) > 0, 


then there is no change in the maximum pressure. The only other possible 
case is 

fo > 9 > (vyM+0—1)/(yM+ 26) 
in which the peak pressure always occurs at ‘burnt’. The change in peak 
pressure is obtained from (29) by writing f = 0, giving 


C RT (1 —fo) )AP,,/ on AP, Ax 


(y—1) fo/{1—Hy—1) M} + 

Z, pi 
(1=fo1—Hy— DM] 
+-(y—1)M{1+6—6,(1—fy) [1 —3(y— 1) M}9— Zp 17/0, (1+8). 


+{1+9—(y— —Fo)}(1+-8)*In 


It is sufficient to test this against exact solutions for concentrate 
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resistances. The agreement is shown in Tables 4 and 5, and is as good as 
the corresponding results for the muzzle velocity (Tables 1 and 2). One 
need not test the accuracy of superposition, which will obviously run 


parallel to the accuracy obtained for the velocity. 


TABLE 4 





CRT,(1—fp) AP m 
P, Ax/7°99 P,, AP, Ax Method 
1°89 First-order theory 
| Exact integration 
5 2°18 sy » 


” ” 





{ccuracy of first-order the ory at 0-58 calibres travel. 
f., 0-9; P, Aa 7-99 ton/in. stops the shot. 


TABLE 5 








AP,,,, exact First-order 
lo inlegralion theory | Error 
1°78 2°35 2°43 
" 0-76 0°03 
7 "33 0-02 
6 I O'14 | —o-03 
Standard bore resistance, Py Aa 7:99 ton/in., at various positions along the bore. 


| am indebted to the Chief Scientist, Ministry of Supply, for permission 
to publish this paper, to Professor N. F. Mott, F.R.S., for his encourage- 
ment of the work, and to Mr. Tadman for arranging the solution of test 
examples on a National machine. 
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FLOW PATTERNS AND THE METHOD OF 
CHARACTERISTICS NEAR A SONIC LINE 
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SUMMARY 

The work described has two aims: (i) to investigate the geometry of the flow 
pattern near a sonic line, and (ii) to develop a method for the numerical integration 
of the equations of motion near a sonic line, where the usual method based on 
characteristics (ref. 1) breaks down. Only isentropic irrotational steady flow with two 
independent variables is considered, and attention is mainly confined to plane flow 

(i) The directions of the Mach lines coincide at a sonic line, and, in general, their 
curvatures are infinite; the Mach lines are approximately semi-cubical parabolas 
(eqn. (21) ), the isoclines approximately parabolas (eqn. (15) ), and the equipotential 
lines approximately cubical parabolas (eqn. (14) ). The state of affairs is different if 
the curvature of the streamline vanishes on the sonic line; the curvatures of the 
Mach lines at the sonic point are then discontinuous but bounded ; other results 
for this special case are set out in § 6. 

(ii) Two methods are described for extending computations near a sonic line. 
In one method double Taylor series for the velocity components are used; in th 
other method expansions in the neighbourhood of a sonic point are used for the 
coordinates of a Mach line through the sonic point in powers of the complement 
of the Mach angle. The methods are developed for plane flow; the necessary 
modifications are indicated for flows with axial symmetry. 


1. Taylor expansions of the velocity components near a sonic point 
TAKE some sonic point as origin of rectangular Cartesian axes with the 
a-axis in the flow direction. We assume that the velocity components are 
analytic functions of x and y; then, denoting the values of quantities at 
the origin by the suffix s and partial derivatives by suffixes x and y, we 
may write 
Rows 2 
U = Us PLU, Hp YUy + > (UU gy + 2eYU gy +Y*Uyy) +... | 
1 " 

| ) 0295 EIT | a295 | 
v= U,+20, bYy +s (% Veet 2LYVpy+Y Yyy) + 


(For clarity the suffices s are omitted from the velocity derivatives. 
A derivative at a general point is distinguished by its full symbol, éu/éx 
for instance.) The equations satisfied by wu and v are 

paiad _— —_ (2) 


Cx Oy 


: >, OU ou , >, OV ‘ 
(u?—a?) — + 2uv — + (v?—a?)— = 0, (2) 
7 7 I a 
Ox ap) ey 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 2 (1949)] 
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where 2a? \Y 1)q? = (y+1)a?. (4) 
With our choice of axes v, = 0, and hence u, = q, = a,, the critical 


velocity. Expressing each term in equations (2) and (3) as power series 
inz and y by use of the equations (1) and (4) and comparing coefficients 


of corresponding powers, we find 


7 u 2 0 
y y ? 
2fy+1 os, ) . 
Je Uny Dice Uyy 2 fy a Uz uy " (9) 
a,\ - J 
“u ul 
\y 1) . etc. (6) 
ty } a. 


We note that the flow can be determined at any point near the origin 
provided only the x derivatives of u and v at the origin are given. 
If h,do is the element of length in the stream direction, «, the curvature 


of the streamline, then 


ul { ob 7 ss 
a fend Ko)s (7) 
cq 
ant u - 8 
- i : (8) 


\ point on the sonic line at which («,), 0 will be called a singular sonic 


point. 


2. The slope and curvature of the sonic line 

To avoid ambiguity in the sign of the curvature of the sonic line we 
first define the positive direction of the sonic line as that obtained from 
the local stream direction by counterclockwise rotation through an angle 





», With 
o6 
0 " 7 When .. » - 0, 
tg CO 
all 
7 7 0 when a. © <2 
i, 00 
| oq . . 
and, when 0, 7 = 0 or w according as x, < 0 ork, > 0. 
t, ¢ CO 


From the equation of the sonic line, g? = a2, with q? expressed in terms 
of and y by equation (1), we obtain the slope of the sonic line at the 
origin 


da 
| © an ——t on tone, (9) 


dx u 
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If h,ds denotes the element of length in the positive direction of the sonic 
line, the curvature of the sonic line at the origin is 


: 6(8+-n) 
“e hes }, 


i . 
a F,4_1 2 
= — {Uy +A, Uy U 


' 9 
— +A, UF Uyy—2 4 
aut u2)i an Tg Uz Uyy A, Uy Uy Urys (10) 


This can also be written, without reference to any particular system of 
coordinates, in the form 


13 
COS"). ‘ ' 
cs {27, a tan 7+ 2x24 tan?y—(y+1)2atan?7—q,,},, (11) 
a,(K,)s : 
OK C O 
where tT =>—,; Veo sales) 
h,0o h, eo\h, 60 


3. The curvature of Mach lines at the sonic line 

The angle between the positive direction of the + Mach line through 
the point (x,y) and the positive x-axis is « = 0-y (ref. 1), and at the 
origin the positive direction of the + Mach line is along the + y-axis. 
Let h,da, hgdB be the elements of length along a Mach line and along 
the normal to a Mach line, respectively. The curvature of a Mach line, 


kK, = €e/h, x is given by the second characteristic equation (1) 
ov Ov ee Vv, Op 
UB log l —< eg B — (0% T 2eg| +—* — == @, 
‘ \ hy x h, Ox P | p hex 


With the help of Bernoulli’s equation this may be written 





. y—1+2 cos 2v OC 
Y 
m : ( 


K ———— ; (12) 
2q?sinveosv ~h,éx’ 
where v = 4n—p. 
Hence, near the origin, on a Mach line through the origin, 
+1 (k,), . 
” _ +! (Ko), 4.O(1). (13) 


At any sonic point the two Mach lines have a common tangent. At all 
non-singular sonic points («,), is infinite; near such points, on both Mach 
lines, x, and («,), have opposite signs. 


4. The behaviour of the equipotential and isocline near an ordinary 
sonic point 

By expanding the velocity potential as a Taylor series and using equa- 

tions (5) and (6) we can show that the equipotential through the origin 

has the form i 

6a x 

y? nn ee (14) 

(y+1)u, Uy, 
near the origin. 
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FLOW 


The equation of the isocline (i.e. the line along which the velocity has 
constant direction), through the origin, is v = 0. Hence, from (1), (5), 
6), the isocline approaches the parabola 


oP a secu (15) 


The relative positions of each Mach line, streamline, isocline, and equi- 
potential through a non-singular sonic point depend on the signs of (q,), 





Fic. 1. (q)s 0, (Kg)s 0. OA and OB are respectively -+- and — Mach lines, 
OS is the sonic line, OE, TOT’ and FOF’ are the streamline, isocline, and 
equipotential respectively. 


and («,), and are easily found with the aid of equations (9), (13), (14), and 
(15). As an example, they are shown in Fig. 1 with (q,), and (x,), both 
positive. 


5. Determination of the flow in the characteristic triangle based on 
a sonic initial line 
When initial data are given on the sonic line the usual method of 
characteristics cannot be used to find the continuation of the flow in the 
supersonic region. Of the two methods here designed for this purpose 
the first employs the double Taylor series (1) to find the flow at chosen 
points in the supersonic region on normals to the sonic line, and the 
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second determines the coordinates of points on both Mach lines through 
sach sonic point corresponding to chosen values of v. Each method enables 
us to construct a new initial line beyond which the continuation of the 
flow can be found by the usual characteristics method. 


5.1. First method of continuation 

At a point in the supersonic region at a given distance, €, along a 
variable curve through an origin on the sonic line, the value of p, and 
hence the value of the Mach line curvature, is least if the curve is the line 
of steepest pressure gradient through the origin. If we replace it, for 
simplicity, by the straight normal to the sonic line through the origin, 
we can determine the velocity components at a point on a new initial line 
by taking a given sonic point as origin and evaluating the series (1) at 
the point (x,y) where 

Xx F y : é 


sin 7, COS 5 


The value of € corresponding to a chosen value of v on the new initial 
line can be found approximately from the equation 





gy 
= = sec*y, 
a 
‘ : . 1 cosy, 2, ‘ ; 
which gives g —_—— “Sy? O(v4). (16) 
ae l Ko)s 


5.2. Second method of continuation 
Referred to Cartesian rectangular axes through a given sonic point with 
the x-axis in the flow direction, the coordinates of any point in the super- 
sonic region are functions of @ and ». At each point on a given Mach 
line of either family a fixed relation exists between 6 and p, given by the 
second characteristic equation; the coordinates of a point and the values 
of flow quantities on a Mach line can therefore be expressed as functions 
of « only. Provided these functions are analytic we may represent the 
coordinates near the origin of either Mach line through the origin by the 
power series a 
x= a,v+a,v?+..., 


y = 6,v+-6b,v?+.... (17) 


When the origin is an ordinary sonic point the slope and curvature ot 


each Mach line are infinite at the origin. Hence a, = 6, = a, = 0. 

The coefficients a5, @4,..., bg, bs... are found in terms of physical constants 
and the boundary conditions at the origin from two identities in v. In the 
first characteristic equation, 

dx 


—_— = —tan(#-+-v), (18) 
dy 
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dx/dy is expressed in terms of v by means of (17) and @ is related to v by 
combining the second characteristic equation with Bernoulli’s equation, 
whence 


23 {1 5—y i 2(y* 28y +61) 


= —_ = = - y*+- — 


AL, 19 
yt1|3" 15(y+1) '  315(y+1)? ' (19) 


The second identity is provided by the strip condition satisfied on a Mach 
line, i.e. the condition that g, x, y are all functions of v only, 
dq eq dx , oq dy 


( ( — —j— (, 
lh i f 


, 20 
ex dv ° * dy dv sates 


‘ ; ; di 
All the terms of equation (20) can be expanded in powers of v: 7 by 
, 


Bernoulli’s equation, q: a< and qa dy by means of (1) and (17). 
ox di dy dv * 

Equating coefficients of corresponding powers of v in equations (18) and 

(20), and solving the resulting equations for dg, d4,..., bg, bg,..., we obtain 

the following expansions in powers of v for the coordinates of the +- Mach 


line through an ordinary sonic point 





| asd Fr 6 
M__ty 5 ue_ _5_ Gy they] ¥* 
2 3 3 “= sels 


er : ae? 2 = 
3(} 1) u, dU, AY ) Uy 6 


- a 3 a4 > 
y+ 10 Uy, , Uz 5d A, Ury v | (21) 
B(y+l)u, Sus B(yt+l1) uw J5 | 


Equation (21) shows that near the origin each Mach line behaves like one 
branch of the semi-cubical parabola 


9a. 9 
2(y+1)u, 2(y+1)(kg)s 


The series (21) diverge when u,/u, is sufficiently large. The first method 
must therefore be used to construct the new initial line when the origin 
lies close to a singular sonic point. When the origin is itself singular the 
Mach lines through it can be represented by the special power series 
described below. The second method is particularly useful for finding 
the continuation of flow in a local supersonic region of potential flow; 
\k,), and uw, are everywhere different from zero on the sonic line bounding 
such a region, and singular points are absent. 
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6. The Mach lines, the equipotential, and the isocline at a singular 
sonic point 
We have seen that when («,), = u,/a, = 0 the curvature of either Mach 
line through the origin may be finite at the origin. In finding the expan- 
sions for the coordinates of such a Mach line we therefore retain the 
coefficients a, and b,. Equating coefficients of corresponding powers of » 
in equations (18) and (20) we obtain two sets of equations for dg,..., 5 


1 Sees 
which now begin with 
2A5 Fb, 
> 2 9 
F 9 2 3 oe 2 
and 24, UA, +A, Uy, b5 = rary hac 
: a 
respectively. Solving these with the aid of (5) we find 
— a, 2a, 
+b, 2a, = ———- or ———+*. 
+)u, (y+)u, 
Hence from equation (12) 
+1)u, _(y+1)u, 
ih xe ge FOI (22) 
x rs 
a, 2a, 


Thus each Mach line has a discontinuity of curvature at the singular 
sonic point and has two branches. From equation (10) the curvature of 
the sonic line itself at the singular sonic point is x, = (y+1)u,/a,. The 


first branch of each Mach line is convex to the subsonic flow with 


Ka)s +x,: the second branch is concave to the subsonic flow with 
a/s a 
(k,)s = +4«,. The equipotential through a singular sonic point can be 
shown to behave, near the origin, like 
24a3 x 
y3 s 


a” a 

(y+1)°uz 
The corresponding isocline has a double point at the singular origin. 
One branch touches the y-axis with curvature }(y+1)u,/a,; the slope of 


the other branch at the origin is —}w,,,/u,,,.. The locusx, = 0 lies between 


yy* 
these two isocline branches, with slope at the origin equal to —w,,,/u 


ry! yy’ 


The complete flow pattern near a singular sonic point at which (7,), 
and (q,), are positive is shown in Fig. 2; it agrees with that found by 
Lighthill (2) in symmetrical channel flow. The + Mach lines represented 
by PQR do not originate on the sonic line; they are the reflections of 
— Mach lines starting from OS. 

It is of interest to determine the point at which the curvature of the 
Mach line PQR changes sign. Craggs (3) has shown that the branches 
OA’, OB’ of the singular Mach lines in Fig. 2 are branch lines. It can be 
shown further that a point where the Mach line curvature vanishes must 
lie on the isobar where the Mach number M = 2/,/(3—y) or on a branch 
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line or a higher singularity of the same type. Since we are considering 
fields of flow where the velocity derivatives are bounded and continuous, 
the Mach line curvature can change sign only where it vanishes (except 


F P 
ie B 


ba . \\ 
~ . 
™ | 
~ 
“i % E 
~ 
™~ 
0 ‘ip 
= 
. ~ 
~ 
1 * ~ 
r\- 
| 
7 
E fl . s 
4 
/ “-g4 R 
/: k 

/] 

/ . 

a 
/ F 
/ / 

S i. 4 
B T F 
A’ 
Fic. 2. (x-)e 0, (t3)s > 95 (Go)s 0. O is the singular sonic point, SOS’ is the sonic line, 


40A’ and BOB’ are singular Mach lines of the + and — family respectively, TOT’ and 

t0t’ are the two isocline branches; HOE’ and FOF’ are respectively the streamline and 
equipotential through O; kOk’ is the locus xg 0. 

at the singular point itself), and it follows that the curvature of each Mach 

line PQR changes sign at Q, the point where it crosses A’OB’. 


7. Completion of a known supersonic field up to its sonic boundary 

Consider a field of supersonic flow which is determined from given 
conditions on an initial line inside the field. As we integrate along Mach 
lines towards the sonic line a stage will be reached at which the Mach 
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line curvature is too large for further application of the usual methods 
of integration by characteristics to be possible. We now give two methods 
of finding the continuation of the flow beyond this stage up to the sonic 
boundary. 


7.1. First method 

Suppose O is a point in the known flow field, beyond which the usual 
method based on characteristics becomes inaccurate. With O as origin 
we take rectangular Cartesian axes with the x-axis towards the sonic line 
and the y-axis in the direction of the isobar at O. If the velocity com- 
ponents are analytic functions of x and y, their values at a point near (0 
on the x-axis are given by 


e l eS 3 
u Ug Uy, x + Si ae Us 
) (23) 
! I 2 
t Vo T Uy U5) Ure UT 
where wu, = (€u/éx),, ete. 


The coefficients u,, u,,,... are easily found if the supersonic field behind 
O has been integrated from isobar to isobar. 
The sonic point on the z-axis is then found by solving the equation 
: E 6 9 
urty? — a? 
for x by successive approximation. 
7.2. Second method 
Using rectangular axes through O, with the z-axis now in the stream 
direction at O, we represent either Mach line near O by expansions 
«= A,h+424¢"+... | 
| , Where ¢ = vp—v, v9 = (v),<9- 
y = b,4+b,¢2+... J 


We express @ and qdq/d¢ as power series in ¢ using the relations 


9 —3 
q Umax f + 7 cost(g—n)| , 
dq 
—tan(v,—d)-+ dé = 0. 


q 


We may then express the first characteristic equation and the strip 
condition as power series in ¢ and by equating coefficients of correspond- 
ing powers obtain two sets of equations giving a,, dy,..., b;, bg,... in terms 
of boundary conditions at O. The sonic point on each Mach line through 
O is then given by 


ce 
L = AyVpt+Ayvo+... | 


9 e 
y = by +b,2+... J 
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FLOW 


8, Application of the methods to flow with axial symmetry 

All the methods developed for plane flow can in principle be extended 
to apply to flow with axial symmetry. Relations between the coefficients 
in the Taylor expansions for uw and v have to be revised; to the expression 
on the left of the governing equation (3) must be added the term 


usin 64+ cos 0, 





ro+asin 0,+y cos 0,’ 


where 7, is the distance of the origin from the axis of symmetry and 6, is 
the angle between the x-axis and the axis of symmetry; otherwise the 
first method of continuation is unaltered. The second method of con- 
tinuation suffers an additional change. The relation obtained from the 
second characteristic equation now involves x and y in addition to 6 and 
. so that, while the coordinates and values of flow quantities on a Mach 
line can still be expressed as functions of » only, the relation between @ 
and p is no longer universal but varies with the boundary conditions at 
the origin. In the development of the second method we have therefore 
to find a third set of coefficients, arising in the expansion of @ in powers 
of v. The three sets of equations now required are found from the first 
and second characteristic equations and the strip condition, all expressed 
as identities in v. Corresponding changes are made in the methods of 
completing a known supersonic field up to its sonic boundary. 

When the origin is taken at the singular sonic point on the axis of 
symmetry the relations between coefficients in the Taylor expansions take 
a particularly simple form and the expansions in powers of v for the 
coordinates on the Mach lines through the sonic point can easily be found. 
The following results are obtained: 


The curvature of the sonic line at the origin = $u,(y+-1)/a,. 


The curvature of the Mach line at the origin is 


HV5+1)(y+1)u,/a, or F-}(V5—1)(y+1)u,/a,. 


/ 


The results should be compared with those found in plane flow. It is seen 
that the flow pattern is qualitatively the same as that in Fig. 2. One 
difference may be noted; the curvature of the Mach line PQR no longer 
changes sign on crossing the line OB’, unless the flow at Q is parallel to 
the axis. 


The author wishes to thank Professor Goldstein and Dr. R. E. Meyer 
for their criticism of this paper. 
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